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ABSTRACT 

We describe a numerical scheme for studying time-dependent, multifluid, 
magnetohydrodynamic shock waves in weakly ionized interstellar clouds and 
cores. Shocks are modeled as propagating perpendicular to the magnetic field 
and consist of a neutral molecular fluid plus a fluid of ions and electrons. The 
scheme is based on operator splitting, wherein time integration of the govern- 
ing equations is split into separate parts. In one part independent homogeneous 
Riemann problems for the two fluids are solved using Godunov's method. In 
the other equations containing the source terms for transfer of mass, momentum, 
and energy between the fluids are integrated using standard numerical techniques. 
We show that, for the frequent case where the thermal pressures of the ions and 
electrons are <C magnetic pressure, the Riemann problems for the neutral and 
ion-electron fluids have a similar mathematical structure which facilitates numer- 
ical coding. Implementation of the scheme is discussed and several benchmark 
tests confirming its accuracy are presented, including (i) MHD wave packets 
ranging over orders of magnitude in length and time scales; (ii) early evolution 
of mulitfluid shocks caused by two colliding clouds; and (iii) a multifluid shock 
with mass transfer between the fluids by cosmic-ray ionization and ion-electron 
recombination, demonstrating the effect of ion mass loading on magnetic precur- 
sors of MHD shocks. An exact solution to a MHD Riemann problem forming the 
basis for an approximate numerical solver used in the homogeneous part of our 
scheme is presented, along with derivations of the analytic benchmark solutions 
and tests showing the convergence of the numerical algorithm. 
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MHD — plasmas — shock waves — waves 

1. Introduction and Motivation 

Shock waves in molecular clouds and star-forming regions are energetic events which can 
process the interstellar gas and dust they entrain. Shock-excited molecular H2O has been 
detected by Herschel (van Dishoeck et al. 2011) and shocked species such as H 2 and CO 
are among the targets of opportunity for SOFIA (Eisloffel et al. 2012). Observations of the 
kinematics and proper motions of protostellar outflows and their associated shocks yield ages 
as young as ~ 10 2 -10 4 yr (e.g., Gueth et al. 1998; Hartigan et al. 2001). Because the time 
scale for these shocks to develop steady flow is typically > 10 4 yr (see time-dependent 
solutions are generally required to model observations of the shock-excited emission. For 
instance Gusdorf et al. (2011) applied approximate time-dependent models to H 2 , SiO, and 
H2O lines in the BHR71 bipolar outflow shocks and concluded that the best fit to the 
observations occurred in models having ages < 2000 yr. However there was a fairly large 
degeneracy in the parameter space which allowed reasonable fits. 

Modeling shocks in weakly-ionized interstellar clouds is a complex task due largely to 
the complex nature of the clouds themselves: molecular clouds are threaded by interstellar 
magnetic fields of order 10-10 3 //G (Crutcher 2004; Rao et al. 2009; Tang et al. 2009) and 
are only weakly ionized, with ion fractional abundances < 10~ 5 (Caselli 2002; Miettinen et 
al. 2011). Because the bulk of the matter is neutral, clouds have a finite electrical conduc- 
tivity; consequently the neutral gas and charged particles collectively constitute a nonideal 
magnetohydrodynamic (MHD) system. Moreover the neutral particles respond indirectly 
to magnetic forces via momentum exchange with the (rare) ions and electrons, so that a 
multifluid treatment of the dynamics is required. Early work on multifluid MHD shocks 
established that the ions can be accelerated ahead of the neutral shock front in a "magnetic 
precursor" (Mullan 1971; Draine 1980). Ion collisional drag in the precursor accelerates, 
compresses, and heats the neutrals, thereby passing information about the disturbance on 
to the neutral gas upstream from the shock. Depending on the shock speed, the time scale 
for collisions between the ions and neutrals, and the rate at which the neutral gas can cool, 
the resulting MHD shock can be a discontinuous flow (J-shock) with a supersonic to sub- 
sonic transition at a shock front, a continuous flow (C-shock) which is supersonic everywhere 
(Draine 1980; Chernoff 1987; Draine Sz McKee 1993), or a continuous flow with a supersonic 
to subsonic transition (C*-shock, Roberge & Draine 1990). 

The majority of theoretical work has been done on steady multifluid MHD shock waves. 
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The chemical abundances and excitation of certain molecular species in shock waves prop- 
agating perpendicular to the ambient magnetic field ("perpendicular shocks") were studied 
by Draine et al. (1983). Steady perpendicular shock models with increasingly elaborate 
physics (including the effects of water emission, extensive chemical networks, etc.) were sub- 
sequently presented by Flower et al. (1985, 1986), Pilipp et al. (1990), Tielens et al. (1994), 
Kaufman & Neufeld (1996a,b), Flower et al. (1996), Schilke et al. (1997), and Guillet et 
al. (2009). Steady MHD shocks propagating obliquely to the magnetic field have also been 
studied (e.g., Wardle & Draine 1987; Caselli et al. 1997). However the applicability of steady 
models is limited when it comes to interpreting observations of shocks in star-forming regions 
for reasons noted above. 

A study of time-dependent, multidimensional C shocks using a two-fluid finite-difference 
formulation, and including the inertia of the ions, was presented by Toth (1994). Time- 
dependent simulations of the formation of C-shocks in models which neglected the inertia 
of the ion-electron fluid were produced by Smith & Mac Low (1997). MacLow & Smith 
(1997) investigated the nonlinear development of Wardle instabilities in three-dimensional 
C-shocks (Wardle 1990) using a two-fluid, time-dependent, finite-difference MHD code, in- 
cluding the inertia of the ions. Ciolek & Roberge (2002) and Ciolek et al. (2004) simulated 
the formation and evolution of one-dimensional perpendicular shocks, accounting for the in- 
ertial effect of charged dust grain fluids (ion inertia was ignored, however). Time-dependent, 
one-dimensional, perpendicular shock models including the inertia of the ions, heating and 
cooling, and an extensive network of reactions for 33 different chemical species, were de- 
scribed by LeSaffre et al. (2004a,b). Multidimensional and multifluid (including charged 
dust grain species, but neglecting ion inertia) evolutionary MHD shock models were dis- 
cussed by Falle (2003) and Van Loo et al. (2009); the same numerical code was also used 
to investigate the effect of upstream density perturbations on dusty C-shocks in Ashmore 
et al. (2010). Finally, a one-dimensional study of the redistribution of mass and magnetic 
flux and potential protostellar core formation due to ambipolar diffusion (the drift of the 
charged particles with respect to the neutrals) in transient C-shocks arising from colliding 
flows was presented by Chen & Ostriker (2012); in their calculation the inertia of the ions 
(and self-gravity of the system) was not included. 

A pseudo-time-dependent (quasi-Lagrangian), one-dimensional MHD shock model has 
been employed in several studies of multifluid shocks (e.g., Chieze et al. 1998; Gusdorf et 
al. 2008; Flower & Pineau des Forets 2010, 2012). In these models, time dependence is 
mimicked by setting the partial derivative d/dt equal to zero in the governing equations 
for conservation of mass, momentum, energy, and magnetic flux (see eqs. [I] - [6] below), 
and then replacing the spatial derivatives d/dx in those equations with a "flow derivative" 
(l/v)d/dt, where v is the fluid velocity of either the neutrals or the ions. The pseudo-time- 
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dependent approach has grown to include increasingly more detailed molecular and chemical 
networks. 

To study MHD turbulence in weakly ionized clouds, some have advocated a multifluid 
MHD method referred to as the "heavy ion approximation" (Li et al. 2006; Oishi & MacLow 
2006; Li et al. 2008). This approximation was introduced to overcome numerical difficulties 
associated with large values of the ion Alfven speed, which are typical in molecular clouds 
(see eq. [M] below) and limit the size of the time steps a numerical code can take to 
produce accurate, stable simulations (via the CFL condition, see eq. [29]) . In the heavy ion 
approximation the masses of the ions are artificially raised to increase the mass density of 
the ion fluid so that smaller, more tractable ion Alfven speeds (and thus, CFL-limited time 
steps) are attained. To keep the collisional drag terms between the neutral and ion fluids 
unchanged, the coupling constants which appear in the ion-neutral frictional forces (e.g., 



see eqs. 13 -[15]) are also adjusted so as to offset the effect the inflation of the ion masses 
has on the forces. However, Tilley & Balsara (2010) showed that, while the collisional force 
terms remain unchanged in the heavy ion approximation, characteristic length scales (which 
involve different combinations of the ion masses and coupling coefficients) associated with 
the dissipation range of MHD turbulence are not calculated correctly using this method. 

Roberge & Ciolek (2007, hereafter RC07) examined the initial phases of MHD shock 
formation in weakly ionized clouds, including the effects of ion inertia. They noted that for 
sufficiently small times the disturbance in the ion-electron fluid is linear, and derived explicit 
analytic solutions for the time evolution of a perpendicular shock. RC07 found that at very 
early times the inertia of the ions determines the formation and propagation of the magnetic 
precursor to a neutral shock. At later times the ions' inertia becomes negligible, their motion 
is force free, and the evolution of the precursor is nearly self-similar. 

The essential mathematical difficulty in modeling time-dependent, multifluid shock 
waves is caused by the transfer of mass, momentum, and energy between the fluids by elastic 
ion-neutral scattering, ionization/recombination, and numerous other atomic and molecular 
processes. In the absence of coupling between the fluids, the multifluid shock problem would 
reduce to independent solutions of Euler's equations for the neutral gas and the equations 
of ideal MHD for the ion-electron fluid. In both cases the governing equations are hyper- 
bolic partial differential equations (PDEs) which can be solved using well-known techniques. 
Many of the algorithms for Euler's equations rely on one's ability to solve the "Riemann 
problem" (e.g., Courant & Friedrichs 1948), wherein the gas initially is in two uniform but 
different states separated by a discontinuity. In particular, Godunov (1959) realized that the 
evolution of gas in two adjacent cells of a finite difference grid is indeed a Riemann problem 
and exploited this fact to produce an algorithm, Godunov's method, whose descendants ac- 
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count for a large subset of all algorithms for gas dynamics (e.g., see Toro 2009 and Pirozzoli 
2010). Riemann solutions are combinations of the characteristic waves of a fluid — sound 
waves, rarefactions, and contact discontinuities in the case of a gas — which are orchestrated 
to satisfy certain matching conditions where the waves intersect (see App.|A]for an example). 
Riemann- Go dunov algorithms also exist for ideal MHD but are complicated by the fact that 
seven characteristic waves exist when the flow an arbitrary angle with respect 

to the magnetic field (see Fig. 2 of Dai & Woodward [1994] or eq. [7] of Torrilhon [2003]). 

When interfluid coupling is included, the only changes to the hyperbolic PDEs for each 
fluid are the addition of certain "source terms" which describe mass, momentum, and energy 
transfer between the fluids (see eq. [l]-[([6]). Toro (2009) noted that problems of this nature 
could be solved by a technique called operator splitting and gave a simple, nonhydrodynam- 
ical example. Operator splitting combines separate solutions of the "homogeneous problem" 
(source terms set to zero) obtained, e.g., with a Riemann- Go dunov alogrithm, and the "in- 
homogeneous problem" (evolution due only to the source terms) obtained with standard 
methods for ordinary differential equations (ODEs). Tilley et al. (2012) have recently pre- 
sented an operator-splitting scheme for multifluid MHD which treats general time-dependent 
flows for a two-fluid system of neutrals and ions traveling at any orientation with respect to 
the magnetic field. Their algorithm retains the inertia of both the neutrals and the ions. Col- 
lisional drag between the two fluids (with a velocity-independent collision rate) is included 
along with an energy equation with source terms for each fluid. The algorithm of Tilley et 
al. (2012) does not account for mass exchange between the neutral and ion-electron fluids, 
presumably because this is not important for the applications of interest to them. Tilley 
et al. (2012) presented various benchmarks tests, including the development of a Wardle 
instability, and found that the scheme performed well. 

In this paper, the first in a series on shocks in molecular clouds, we present a split- 
operator method similar to the algorithm of Tilley et al. (2012), but tailored to the one- 
dimensional geometry appropriate for shock waves. This allows us to exploit the unique 
property of perpendicular shocks in weakly ionized plasmas: there exists an exact solu- 
tion to the MHD Riemann problem, which is fundamental to our algorithm. While the 
assumption of perpendicular shocks is obviously a special case, it does treat the fundamental 
mathematical problem of coupled hyperbolic PDEs including mass, momentum, and energy 
transfer. Deferring the case of oblique shocks to future work also makes strategic sense be- 
cause modeling perpendicular shocks has a particular advantage: for likely cloud and shock 
parameters, the thermal pressure of the ion-electron fluid can be neglected compared to 



magnetic pressure (§ 2.1). In this regime the charged fluid acoustic modes play no part in 
the MHD Riemann problem and the number of characteristic waves in the ion-electron fluid 
reduces from seven to just three — the same number of waves that occur in the gas dynamic 
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Riemann problem. As a result, it turns out that the Riemann problem for the charged fluid 
can be solved exactly (App.[A]). From this exact solution an approximate but accurate solver 
can be constructed for use in the numerical solution of the MHD Riemann problem in the 
split-operator method. In our algorithm we also include the effect of mass transfer between 
the neutral and charged fluids, which can have profound effects on multifluid shocks (Flower 
et al. 1985). 

The plan of our paper is as follows: in § 2 we present a formulation of the model, 
including the governing equations and assumptions, the respective Riemann problems for 
the neutral gas and ion-electron fluid, and a description of how the split-operator scheme 
is implemented to evolve model shocks in time. We assess the accuracy of our numerical 
solutions with various benchmark tests in § 3. Our results are summarized in § 4. In the 
appendices we present the exact solution to the MHD Riemann problem for perpendicular 
flows (App. |Aj), an approximate solver based on this solution (App. [B]), a derivation of the 
analytic solutions used as benchmark tests (App.lC]), and tests confirming the order of spatial 



and temporal convergence of our numerical algorithm (App. D). 



2. Formulation 

Our area of interest is shocks and related flows in interstellar molecular clouds. These 
clouds contain neutral particles with number density n n and particle mass m n plus a weakly 
ionized plasma of singly-charged ions and electrons having number densities n> x and n e , and 
particle masses m ; and m e , respectively. Since we are not concerned with the problem of 
gravitational collapse and deal with systems typically having length scales much smaller than 
the Jeans length, self-gravity of the gas is ignored. We adopt a cartesian coordinate system 
(x,y,z) and restrict our attention to models in which all of the physical variables are functions 
of time and the x-coordinate only. We assume that there is a magnetic field B = B(x,t)z, 
and that each of the fluids has a velocity v a = v a (x, t)x, with a = n, i, e. Thus we consider 
perpendicular shocks only; however we note that the split-operator method described in this 
paper can be extended to other geometries. 

In the molecular cloud and core environments we study, the magnetic field strengths 
and neutral gas densities are such that the ion and electron fluids each have Hall parameters 
(= charged particle gyrofrequency times the mean collision time with the neutral gas) ^> 1. 
This means that the ions and electrons gyrate about a magnetic field line many times before 
suffering a collision with a neutral particle, and can therefore be considered to be attached 
to the field. The magnetic field is thus "frozen into" the charged fluid, and the electrons and 
ions move together with v c ~ V\. 
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Finally, we ignore the effects of interstellar dust grains. Grains can become charged 
and numerical simulations have shown that under certain conditions they can affect the 
evolution of MHD shocks (Ciolek & Roberge 2002; Ciolek et al. 2004; Van Loo et al. 2009); 
the influence of dust grains in MHD flows using a split-operator method will be considered 
at a later time. 



2.1. Governing Equations and Assumptions 

Comprehensive derivations of the multifluid MHD equations for astrophysical flows are 
presented in Draine (1986) and Mouschovias (1987). For the geometry adopted here they 
are 

dp n d{p n v n ) 

op a i 

-^ + — ([E n + P n ]v n ) = F Q v n --S n v n 2 + G n -A n (3) 
dpi d(piVi) _ 

dt dx ~ { ' 

+ - -F. (5) 



dt dx \ 8n t 

dB djBvQ 

-dl + ^x~ = ° ■ (6) 

Equations @-@ 

express the conservation of mass, momentum, and energy for the neutral 
gas. Equations (|4])-([5]) express mass and momentum conservation for the ion-electron fluid 
and ^ is the induction equation. We also impose macroscopic charge neutrality, 

e{rii - n e ) = . (7) 

To the set above one should generally add separate energy equations for the ions and 
electrons, which are needed to calculate the ion and electron temperatures, T\ and T e . How- 
ever if elastic collisions dominate the transfer of energy between the ions and neutrals, as is 
usually the case, then 

T i « T * + Sr to - v «? ( 8 ) 

3fc B 

(Chernoff 1987), an approximation we adopt. There is no analogous approximation for 
T c because energy exchange between the electrons and neutrals is dominated by complex 
inelastic processes such as electron impact excitation and ionization. We temporarily omit 
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the electron energy equation and estimate T e when it is needed in one of our benchmark 



calculations (see £3.3) 



The quantities p n (= m n n n ) and p\ (= mini) are the mass densities of the neutral gas 
and ion-electron fluid (we neglect the electrons' contribution to p\). The neutral gas has 
thermal pressure 

Pnk B T n 



n 



(9) 



where T n is the neutral temperature and fee is the Boltzmann constant. Its energy density is 

E n = ^p n V n 2 + P n e n , (10) 

where e n is the internal energy per unit mass. In general it is necessary to find e n by a kinetic 
calculation of the level populations of rotationally and vibrationally excited H 2 , CO, H 2 0, 
etc. (e.g., Flower et al. 2003). Here we temporarily forego this complication and use the 
internal energy for an ideal gas in thermodynamic equilibrium, 

P 

(11) 



" Pn(7"l) ' 

with 7 = 5/3 in the calculations presented below. However it is important to note that the 
split-operator method does not preclude a kinetic calculation of e n or nonideal equations of 
state for e n . 

In the ion-electron momentum equation (|5j) we have dropped the thermal pressure force 
because it is typically much smaller than the magnetic pressure force. The ratio of the 
ion+electron to magnetic pressure is 

-MU + T.) 1.39 x 10-» ( (—^) (At) i 7 ^) , (12) 



(B 2 /8n) \ B J V10 4 cm" 3 / V10- 7 / V 10 3 K 

where X\ = n\/n n is the fractional ionization and we have normalized quantities to represen- 
tative values for a shock wave. We conclude that thermal pressure of the ions and electrons 
is indeed negligible for the conditions of interest here. 

The source terms S n , F n , and G n — A n are the net rates per unit volume at which mass, 
momentum, and thermal energy are added to the neutral gas, where A n is rate of radiative 
cooling and G n is the net rate of heating by all other processes. In £j3] we describe bench- 
mark calculations with various assumptions about the source terms. In all cases momentum 
exchange is dominated by elastic ion-neutral scattering with 

F n W F n ,cl = - — (fn - Vi) . (13) 



The time scale for drag to accelerate the ion-electron fluid is 



(mi + m n ) ( _ 



x>\Vn ~ ^i| 
(™)in 



n 2 



-1/2 



(14) 



where (crw) in is the Langevin collision rate (Giousmousis & Stevenson 1958; Flower 2000) 
and cjgeo is the geometric cross section. It follows from Newton's 3rd Law that the time to 
accelerate the neutral gas is 

Tni = (pJpi)T m > T in . (15) 

The ion-neutral and neutral-ion drag times, T in and T ni respectively, are fundamental time 
scales for the flow. 



2.2. The Spit Operator Method 



The system of governing PDEs (JT|-((6]) has the conservative form 



dU dJ= 

dt dx 



(16) 



where U = [U n , £/;] T is the array of conserved dependent variables, T = [J- n , J-\\ 1 is the 
corresponding array of fluxes, and S = [S n , «Si] T is the array of source terms, with 



t . 



U n = [p n , PnVn , E n ] T , 

Ui = [pi , PiV U B] T , 

= [Pn^n , Pr^n + Pn , + P n )v n ] T 

T\ = [piVi , pM 2 + B 2 /8n, Bvi] T , 

S n = S n , F n , F n W n - -S> n 2 + G n - An 

5i = [S a , - F n , 0] T . 



(17) 
(18) 
(19) 
(20) 

(21) 

(22) 



There are several different ways to solve system ( 16 ). In this paper we apply a scheme known 
as operator splitting, an authoritative discussion of which can be found in Toro (2009; see 
Ch. 15). The basic idea is to split the full problem into two parts, called the homogeneous 
and inhomogeneous sub-problems. In the former one solves the hyperbolic system 
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on a discretized computational domain using, say, a Riemann-Godunov (RG) solver (see 



£2.3). In the latter one solves 

using, say, a Runge-Kutta algorithm for ODEs. 

The accuracy of the resulting solution depends on the manner in which the sub-problems 
are coordinated. For example, U(x,t) can be advanced from time t to t + At with first-order 
accuracy in At by setting 

U(x, t + At) = Tj^T^ [U(x, t)} , (25) 



where the operator Trq ^ advances the solution from time t to t + At by solving (23). The 

result is then fed to Tj ^\ which advances the solution from t to t + At by solving (24). 
Second-order accuracy can be attained by using what is sometimes referred to as "Strang 
splitting" (Strang 1968; Toro 2009), 

U(x, t + At) = T (At/2) r (At) r (At/2) t)] ^ (2g) 

where the inhomogeneous sub-problem is advanced by a "predictor" half step At/2, the 
homogeneous sub-problem is advanced by a full step At, and the inhomogeneous sub-problem 



is integrated again over another half step. We use Strang splitting in our algorithm (see { 2.3 ). 



In the operator-splitting method, the homogeneous subproblem (23) can be further 
separated into independent Riemann-Godunov problems for each fluid: 

" ° • 

This means that the neutral gas and ion-electron fluid do not directly influence one another 
in either of their separate RG problems. The two fluids are therefore dynamically uncoupled 
during this stage of the calculation. Operationally, this has the advantage that existing well- 
developed numerical techniques (such as approximate Riemann solvers for gas dynamics) can 



be applied readily to RG problem (27) for the neutral gas. In Appendix B we describe an ac 



curate MHD Riemann solver for RG problem (28) for the ion-electron fluid. Thus there is an 



overall symmetry to the solution of the two separate Riemann-Godunov problems when the 
operator splitting scheme is used. Exploiting this symmetry greatly enhances the efficiency, 
and simplifies the writing, of a numerical code designed to model multifluid shock waves. 



Another advantage of our scheme is that by having RG problem ( 28 ) as one of the governing 
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equations, the inherent MHD hyperbolic (i.e., wave) structure is unambiguously built into 
the time evolution. That is, the inertia of the ion-electron fluid is always included. Of course 
the neutral gas and ion-electron fluid do influence one another. Dynamical recoupling of the 
fluids occurs during the source integration stage (24) of the split-operator method. As we 
show in our benchmark tests (§ 3), the interaction between the fluids is indeed accurately 
accounted for during this stage of the calculation. 



2.3. Outline of the Algorithm 

Models are calculated on a spatial domain consisting of a set of N fixed mesh points 
{xj,j = 1, • • • , N} with uniform spacing Ax = Xj — Xj-±. Variables are calculated at Xj. 
Located midway between cell j (centered about Xj) and cell j + 1 (centered about is 
the cell face x i = Xj + Ax/2. 

To advance a model in time, a stable numerical time step At has to be chosen. One 
possibility is to impose the Courant-Friedrichs-Lewy (CFL, e.g., Toro 2009) condition at 
each mesh point j, 

v Ax 

A*cfl, = t , (29) 

^maxj 

where A maxJ ' is the maximum (for all wave modes including both fluids) wave speed at grid 
point Xj and v is a number such that < v < 1. In our application it is typically the case 
that the MHD wave speeds far exceed those of the neutral gas. Therefore we set 

Amaxj = (H + V ims )j , (30) 

where V^ ms is the ion magnetosound speed. When ion-electron pressure is neglected, V ims = 
Via, where the ion Alfven speed 

ViA = — 7=p= — 6.90 xtff SU^ff^ (^1) " 2 km s- (31) 



y/Anpi \50 /iG J \ mi J \ X\ J \ n n 

(m p is the proton mass). However the source terms <S n and «S; also contain important time 
scales, including the collisional drag times, r in and r ni , as well as scales related, e.g., to the 
heating and cooling of the gas. Let At$j be some fraction of the smallest time scale associated 
with the source terms at Xj. To have a stable numerical time integration throughout the 
entire computational mesh, it is then necessary that the time step be such that 

At<mm[At Sj ,At CFLj ] . (32) 



Once a stable step size has been determined, integration proceeds as given by equation (26) 
in the fashion described in the subsections below. 
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As discussed early on by Paleologou & Mouschovias (1983, see their App. A), large 
disparities exist between the magnitudes of flow time scales and collisional time scales in the 
multifluid MHD equations rendering the governing equations mathematically "stiff". This 
means that At will often have to be less than the smallest physical time scale in a system 
containing other time scales which are much greater in size. In this circumstance a large 
number of small time steps will then be required to follow a state that is evolving on a much 
longer natural time scale. An extreme example occurs in the third benchmark test presented 



in section § |3.1.3| which is followed in its development until several neutral-ion collision times 
r n i have elapsed, up to a time ~ 2.6 x 10 5 yr (~ 10 to 100 times greater than the ages of 
shocks in star-forming regions that we intend to study — see §[!]). To stably integrate that 
model to that time using an explicit method, At had to be kept below the ion-neutral time 
scale r in « 1.2 x 10~ 2 yr; At = 0.4r in was actually used for that model. Thus, ~ 8 x 10 7 time 
steps were taken to reach completion. While this is indeed a large number of computational 
steps, it is not prohibitively so with modern computational methods. 



2.3.1. Step 1: First Source Integration 



Given some initial values U(x,t), system (24) is advanced in time by At/2. This can be 
carried out using well-known and tested numerical integration methods such as Runge-Kutta 
schemes, or explicit multistep schemes, etc. (e.g., see Ch. 6 of Atkinson 1989; or Ch. 16 of 
Press et al. 1996). For the benchmark calculations presented in § 3, we used a second-order 
Runge-Kutta integrator. The result of Step 1 is an updated state vector U 1 = \U l a , Ul] T 
which becomes the input for Step 2. 



2.3.2. Step 2: Riemann-Godunov Integration 



In Step 2 the dynamically decoupled RG problems (27) and (28) are advanced a full 



time step. The solution for the conserved dependent variables at the end of this integration 
step, U 2 = \U\ , U\Y i is given by 



Ul{x 3 ) 



-—(^Fr(x. i) — T&(x. 

Ax V PV 3+2 J 



(33) 



with 



n,i. 



The variables are known at the grid points, Xj, but the fluxes in equation (33) are 



required at the cell faces. Second-order accuracy in space and time is attained by first 
"reconstructing" the data and then interpolating the variables from the grid points to the cell 
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faces using total variation diminishing (TVD) slope limiter methods such as those employed 
in the MUSCL-Hancock scheme (Toro 2009). For any dependent variable q in the cell 
centered about Xj, reconstruction and interpolation to the cell faces are carried out by 
setting 

1— 1 — 

q(xj) + -Aj , q) = q(xj) - -Aj , 



qf 



(34) 



where R refers to the right face of the cell j (i.e., x. i — S, S — > 0) and L refers to its left 

3+ 2 

face (i.e., x . i + S, S — > 0). Aj is the "limited average slope value" of q. There are many 

3 ~2 

different kinds of slope limiter. For example, the MINBEE slope limiter is 



A, 



max 



mm 



0,min(A i,A / 

3 2 J+ 2 

0,max(A i,A i; 

3 2 J+ 2 



if A i > 0, 

i + 2 

if A i < 0, 

3 + 2 



A , , I = - Qfa) , A 1 = g(^) - 

■ ? " l "2 J 2 



(35) 



(36) 



(Toro 2009), and a generalized version of the monotonic slope limiter of van Leer (1979) is 

1 



(7j mm 



e a 



IA 



i+7 



+ A. 



e a. 



1 < e < 2 , 



1 if A 1 > , a 1 > , 

J ~2 3+ 2 

-1 if A 1 < , A 1 < , 

3 ~2 3+ 2 

otherwise. 



(37) 
(38) 



We normally use the van Leer limiter but both are easy to implement and yield good re- 
sults. We also impose the positivity conditions of Waagan (2009) to the slope limiting and 
interpolation procedure for the neutral gas, which ensures that p n , P n , T n and e n are always 
> at all Xj, even for very hypersonic flows. 

Another TVD-limiter which can be used for the data reconstruction of the charged fluid 
and magnetic field variables is the one derived by Cada & Torrilhon (2009) It has the form: 



qf 



q(x) + ^(T^-)A j+ i , <% = q(x)-±(KKj 1 )* j _i 



where 



and 



max 



Kj 



0,mm , max 



A. 1 

3 + 2 



-0.5ft, min I 211, ^4^, 1-6 



(39) 



(40) 



(41) 
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for flows with smooth minima and maxima, we have found that this reconstruction method 
and limiter function produces smaller relative errors, including the regions about the extrema, 



when compared to that which results from (34) with the limiters (35) and (37). The Cada 



& Torrilhon (2009) scheme is used in the eigenmode convergence test models presented in 
Appendix [Dj 



Integrating equations (27) and (28) over a virtual (predictor) time increment At/2 



further improves the left and right cell face values in each cell: 

At 
2Ax 



< = Uh - ^ (^(Ul) - J=,{U\ 3 )) (42b) 



(Toro 2009). 



=R 



Data reconstruction and interpolation yields variable pairs {U •, U g at x . i and 

\U R ;i, U a •} at a; i, which can be used to calculate the fluxes J-r(x i) and J-b(x i) 
p,j p,jj j-j^ M i+9 H j-~ 



2 ■> 2 



needed in equation (33). Calculation of the fluxes is accomplished by solving the Riemann 
problem at the cell faces using the method of Godunov (1959; for an extremely comprehensive 
discussion see Toro 2009). Solution of the RG problems at each cell interface is carried 
out using approximate Riemann solvers, one for the neutral gas and another for the ion- 
electron fluid. The approximate MHD Riemann solver used for the charged fluid is derived 
in Appendix |B} For the neutral gas we use the approximate gas dynamic Riemann solver 
described in § 5 of Almgren et al. (2010; reportedly based on unpublished work by P. Colella, 
1997). The latter solver is similar to the approximate Riemann solver presented in Toro 
(2009, Ch. 9), but has been extended to RG problems that also include nonideal gases (e.g., 
Colella & Glaz 1985). Using the solver of Almgren et al. allows the modeling of fluids 
having nonideal equations of state, or those that are not in LTE. Our approximate gas 
dynamic Riemann solver is very similar in structure to our approximate MHD Riemann 
solver. (Hence the general symmetry of the problem when using the split-operator method, 
as noted above.) 

The approximate gas dynamic and MHD Riemann solvers are used to obtain the (Go- 
dunov) fluxes at the faces of each cell, J-p{x. i) and J-p{x. i). The details of how this is 

■7+2 J ~2 

done for the ion-electron fluid MHD problem are presented in Appendix [Bj and the details 
for the neutral gas problem follow by analogy (or, see § 5 of Almgren et al. 2010). These val- 

1 



U 2 g are the input for Step 3. 



ues are then inserted into equation (33), thereby giving U 2 3 and U 2 . The updated variables 
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2.3.3. Step 3: Final Source Integration 

Another source integration of size Ai/2 is performed as in Step 1. This yields the 
last update of the physical variables, U 3 = [U 3 , U 3 ] T . The split-operator integration 
scheme is thus completed with the second-order accurate (in space and time) final result 
U(x,t + At) = U 3 . A schematic of our algorithm is presented in Figure [l] Because of 
the modular and general nature of the algorithm we have described, altering or tailoring a 
multifluid MHD split-operator code to one's specific needs, if necessary, would not be difficult. 
For instance, if one wishes to use a different method for interpolation of the variables, such 
as a piecewise parabolic method (PPM; Colella & Woodward 1984), or a weighted essentially 
nonoscillatory method (WENO; Liu, Osher, & Chan 1994), instead of the MUSCL-Hancock 
TVD scheme we describe, the alternate method can be substituted in the indicated data 
reconstruction substeps for both the neutral gas and ion-electron fluids. 

3. Benchmark Tests 

In this section we present some numerical calculations which establish the validity of 
the split-operator method for multifluid MHD shocks and other flows in weakly ionized 
interstellar clouds. The accuracy of the algorithm is tested by comparing our numerical 
results to accurate analytic solutions. Data regarding the convergence and second-order 
scaling of our algorithm are presented in Appendix [Dj In all of the test models the neutrals 
are assumed to be an ideal gas of H 2 molecules having a ratio of specific heats 7 = 5/3, 
and the ion mass m; is set to a generic value of 25m p , which could represent ions such as 
nCO + and Na + . For momentum transfer by elastic ion- neutral collisions we use a Langevin 
collision rate (aw)i Q = 1.7 x 1CT 9 cm 3 s _1 (McDaniel & Mason 1973) and geometric cross 
section <r geo = 7r(r H2 + r ; ) 2 = 2.86 x 1CT 15 cm 2 (r H2 and r; are the molecular and ionic radii, 
respectively). 

3.1. Linear Wave Packets 

In this test we set the source terms S n , G n , and A n to zero in equations ([l]) - ^6|) 
and assume that momentum transfer is entirely due to elastic ion-neutral scattering. We 
initiate small-amplitude disturbances in the plasma and follow the subsequent evolution 
of the charged and neutral fluids. As discussed in Ciolek & Roberge (2002, § 2.3.2) and 
Mouschovias et al. (2011, § 3.2), there exist two important length scales relevant to MHD 
waves propagating perpendicular to the magnetic field. The first length scale is the ion 




Step 3: 
Source Integration 



u 



J 

U(x,t+At)=U 3 



Fig. 1. — Schematic diagram of our split-operator algorithm to simulate multifluid MHD 
shocks and related flows in interstellar clouds and cores. 
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magnetosound wave upper cutoff L ims = 47rVi ms r in . Ion magnetosound waves with wave 
speed V ims = V iA can propagate when their wavelength is less than L ims . The second length 
scale is the neutral magnetosound wave lower cutoff length scale L nias = ttVoa T n i/V nmB . 
Magnetosound waves in the neutral fluid travel at the neutral magnetosound speed, 

V WB = (C n 2 + KA 2 ) V2 , (43) 



o D \ 1/2 / n\l/2 / \ 1/2 / rji \ 1/2 /o X 1/2 

dP n \ (lPn\ nnm / 7 1 / T \ f 2 m p \ _ x 



where 

C ^{^),r^) ^HvaJ [Wk) l^fj ^ (44 » 

is the adiabatic sound speed (the derivative is taken at constant neutral entropy s n ), and 

B ( B \ (2 m p \ 1/2 /10 4 cm- 3 \ 1/2 , . , x 

V nA = —== = 0.771 — — p - km s^ 1 (45) 



y/4:np n \50 /j,G J \ m n J \ n 

is the neutral Alfven speed. Neutral magnetosound waves can propagate for wavelengths 
greater than L nms . At wavelengths between L ims and L nms there is no MHD wave propagation; 
any mode excited on these scales produces ambipolar diffusion of the ions, electrons, and 
magnetic field through the neutral gas. 

To selectively excite different wave modes we impose Gaussian initial perturbations in 
the ion density and magnetic field of the form 

„2~ 



B(x,0) = B 



— X 

1 + A p exp ( — -2 



(46) 



/ x B(x,0) 
ni(x,0) = n i0 ^-^- . (47) 

-DO 

The charged fluid is initially taken to be stationary with Vi(x, 0) = 0. The neutral fluid is 
initially uniform and at rest, with n Q (x, 0) = n n o, P n (x, 0) = P n o, and v n (x, 0) = 0. Quantities 
with a "0" subscript are constant. A p is the dimensionless amplitude of the perturbation, 
and Lq is its width; wave modes excited by this perturbation will have wavelengths clustered 
about L G . 

Below we present three benchmark tests having different packet widths, Lq, but the 
same perturbation amplitude A p = 0.01. All three tests have n n o = 2 x 10 4 cm -3 and 
^io — ^io/ n no — 3.16 x 10 _8 . All particle species are assigned the same reasonable but 
somewhat arbitrary temperature: T n0 = T i0 = T c0 = 10 K. The unperturbed magnetic field 
strength is B = 50 pG, yielding an ion magnetosound speed V ims = Via = 868 km s _1 , 
and a neutral Alfven speed V n A = 0.545 km s _1 . The sound speed C n = 0.262 km s _1 and 
the neutral magnetosound speed V nms = 0.605 km s" 1 . It follows that r in = 1.26 x 10~ 2 yr, 
r ni = 3.19 x 10 4 yr, L ims = 4.33 x 10 14 cm, and L nms = 1.55 x 10 17 cm. 
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3.1.1. Wave Packet 1: L G = 9.35 x 10 11 cm 



Because this wave packet has Lq <C L ims , we expect that the initial perturbation will 
generate two traveling wave packets in the ions and the magnetic field, one traveling leftward 
from the origin and the other rightward. Although acoustic waves in the neutrals can exist on 
these scales (see, e.g., Fig. 1 of Ciolek et al. 2004 or Fig. 6 of Mouschovias et al. 2011), they 
will not be excited by this perturbation: the neutrals are initially uniform and motionless 
and the neutral-ion drag time, r m , is much longer than the decay time scale T& ec for the 
waves (see below). We expect that in this model the ion and magnetic field wave packets 
will propagate through a stationary neutral fluid. 

Given the assumption of stationary neutrals, analytic expressions for the small-amplitude 
wave packets in the ion-electron fluid, which are expected to occur in this model, can be 
derived from the governing equations ([l])-((6]) in the limit where the wavelengths (~ Lq) of 
the Fourier components comprising the packets are <C -^ims- For the initial conditions (46), 



(47), and t>i(x,0) = 0, the solutions are 

i t 



Bi mB (x,t) = B 



1 + ~2 6XP 



27i, 



exp 



[x 



' ims " J 



+ 



ttA p L g 



^i.ims v^J t") 



A 



8 Vj ms Tij; 
-Sims {X-i t) 

V ims exp 



exp 



27i, 



2Ti, 



exp 



erf 



X -\- V; ms t 



J G 



+ exp 
— erf 



L G 
L G 



— exp 



(x + V ims tf 



,(48a) 
(48b) 
(48c) 



(see Appendix C.l). As expected, the solution consists of two traveling Gaussian packets, 
with propagation velocities ±V[ ms and amplitudes A p /2. Momentum exchange (i.e., friction) 
with the fixed neutrals causes the pulses to decay exponentially on a characteristic time 
scale r dec = 2r in . The decay time exceeds r in by a factor of 2 because of the equipartition 
of magnetic and kinetic energy in the ion magnetosound waves which make up the wave 
packets. Elastic collisions with the neutrals directly reduce the kinetic energy of the ion- 
electron fluid; the energy stored in the magnetic field is reduced indirectly as it is gradually 
converted into ion-electron kinetic energy. 



Figure [2] shows the results for this model at t — 5.42 x 10 3 yr using our fully nonlinear 
split-operator method code on a mesh with 5000 grid points. Also shown in each panel 



are the analytic wave solutions (48a)-(48c) evaluated at the same time. The split-operator 



code results are in excellent agreement with the analytic solutions, with the relative error in 
the ion density having a maximum value of 2.5 x 10~ 4 . The split-operator scheme we have 
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Fig. 2. — Model results for the initial perturbation (46), (47), with Lq = 9.35 x 10 11 cm and 
A p = 0.01, shown at time t = 5.42 x 10~ 3 yr. Squares indicate initial values for the neutral 
gas in each panel. Initial ion and magnetic field curves are crosses. Neutrals at time t are 
the solid lines, and the ions are circles. For clarity some data points have been omitted. Also 
shown (dash-dot lines) are the analytic wave solutions (48a)-(48c). Top: number density of 
the neutrals and ions, normalized to their values in the unperturbed state. Middle: fluid 
velocities. Bottom: magnetic field, normalized to its unperturbed value. 



-20- 



employed has accurately reproduced all of the fundamental aspects of the physical evolution 
for this model. This includes the propagation of the magnetosound pulses, the dependence 
of their speed on the magnetic field strength and inertia of the ions (inherent to the solution 
from the RG step of the integration), as well as the decay of the two pulse peaks caused 
by ion- neutral friction (which can only result from the source integration steps). Although 
the ions and neutrals are dynamically decoupled during the RG stage of the integration, 
recoupling and realistic interaction of the two fluids are faithfully reproduced during the 
source integration stages of the algorithm. 



3.1.2. Wave Packet 2: L G = 3.74 x 10 15 cm 



This wave packet has L ims C Lq < L nms . As a result, there will be no propagating 
MHD waves in this case. Instead, the initial perturbation will give rise to an ambipolar 
diffusion mode, in which the ions, electrons, and magnetic field diffuse outward from the 
origin through a neutral fluid which is still effectively stationary on the length and time 
scales characterizing this particular model. The motion of the ions is essentially inertialess 
or force-free; that is, the driving magnetic pressure force is almost exactly balanced by the 
retarding ion-neutral friction ([5]). The linearized ambipolar diffusion mode is found to have 
the analytic solutions 
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(a derivation is given in Appendix C.2). The quantity 

2, 



V ad EE V, 



lms 1 in 



(49d) 



(49a) 
(49b) 



(49c) 



is the diffusion coefficient of the ions and magnetic field through the neutrals. The pulse 
diffuses on the ambipolar diffusion time scale r dec = r ad = L G 2 /4P ad = 37.0 yr 
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Figure [3] shows the results of the split-operator method calculation, in the same format 
as in Figure [2j The output time is t — 174 yr. Also shown in each panel are the analytic 
ambipolar diffusion solutions (49a)-(49c) at that same time (dash-dot curves). The split- 
operator numerical results are seen to be in superb agreement with the analytic solutions, 
having a maximum relative error in the ion density equal to 1.7 x 10~ 5 . The split-operator 
algorithm has a step (the RG step) in which both the neutral and ion-electron fluids are 
described by hyperbolic PDEs, i.e., PDEs which describe wave propagation. In contrast, the 
solution in Fig. [3] exhibits diffusion, which is described by parabolic PDEs (Courant Sz Hilbert 
1953; Dennery & Kryzwicki 1996). The transition from hyperbolic to parabolic behavior as 
the packet width increases is made possible by the source term for momentum exchange 
between the ions and neutrals, which comes into play during the source integration steps 
(Fig. [TJ. The combination of source- and RG- integration steps preserves the underlying 
physics, which dictates that no propagating MHD waves should exist for wave packets with 
dimensions ~ L G if L ims <L G < -^nms- 



3.1.3. Wave Packet 3: L a = 2.24 x 10 17 cm 



This wave packet has L G > L nms , which means it will excite low-frequency (< l/r ni ) 
neutral magnetosound waves over sufficiently large length and time scales. In these waves, 
magnetic forces on the ions are transmitted to the neutrals by collisional drag, allowing 
the neutrals, ions, and magnetic field to participate in collective wave behavior. Because 
Pn ^ Pi ; the overall inertia of the collective motion is due to the neutral fluid. The wave 
speed is therefore the neutral magnetosound speed V^ ms (eq. [43]), with the bulk neutral 
fluid acting as if it is responding directly to magnetic forces on these scales. 

For wavelengths > L nms the ions can still be described as being in force-free motion 
because L nras ^> L ims . Using this fact, we find that the Fourier modes for perturbations with 
L G ^> L nms yield the analytic solutions 
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Fig. 3. — As in Fig. [2] but for Lq = 3.74 x 10 15 cm, and time t = 174 yr. Also shown in each 
panel (dash-dot curves) are the analytic ambipolar diffusion solutions (49a)-(49c) at time t. 
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where 



g x = [1 + (C n /K m s) 2 ] 1 [l+(2P ad t/L G 2 )] V2 
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(51) 
(52) 



(for a derivation, see Appendix C.3). Appearing in these expressions is the neutral thermal 
diffusion coefficient 



V. 



th 



Km' 



V 



ad 



1 + (C n /V nA ) 



(53) 



The last equality follows from equations (49d), (15), (31), (45), and (43). The analytic 
solutions (50a)-(50d) show that there will be Gaussian wave packets propagating with ve- 
locities ±Vn ms , verifying that it is the inertia of the neutral fluid which determines the rate 
at which the packets propagate. The solutions also reveal that the neutral magnetosound 
pulses decay by ambipolar diffusion on a time scale r doc = 2r ad = L G 2 /2D ad = 2.66 x 10 5 yr; 
the decay time is a factor of 2 longer than the diffusion time scale because, once again, there 
is equipartition of magnetic and kinetic energy. 

Figure [4] shows the results of the split-operator algorithm at time t = 2.55 x 10 5 yr 
calculated on a mesh with 2000 points. The figure shows that there are indeed two oppositely 
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Fig. 4. — As in Fig. |] but for L G = 2.24 x 10 17 cm and time t = 2.55 x 10 5 yr. Top: 
neutral and ion densities. Also shown are the analytic solutions (50b) and (50d) for the 
ions (dash-dot) and the neutrals (x). Middle: neutral and ion velocities. Included are the 
solutions (50c) and (50e) for the neutrals (x) and the ions (dash-dot). Bottom: magnetic 
field. The solution (50a) is also displayed (dash-dot). 
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directed wave packets emerging from the initial disturbance about the origin. There is a 
decrease in the neutral gas density from the central region as the matter once located there 
is set into motion and becomes "scooped up" by collisions with the ions, which are being 
driven outward by magnetic field pressure gradients. The transported neutral gas piles up 
just outside the depleted central region and increases the density in the two adjacent pulses 
moving away from the origin. Note that the ion fluid velocity slightly leads that of the neutral 
fluid. This effect is caused by the inertia of the neutrals, which delays the acceleration of 
the neutrals from rest up to the same velocity as the ions. Figure [4] also shows the analytic 



solutions (50a)-(50d) at the same output time. The numerical results are in very good 
agreement with the analytic solutions, with a relative neutral density difference < 3.1 x 1CT 4 
and magnetic field relative difference < 2.2 x 10 



-l 



Because the initial state of the neutral fluid is static and completely lacks density or 
thermal pressure gradients, the RG integration step for the neutrals does not initiate the wave 
motion in the neutrals in this example. Rather it is the collisional drag of the magnetically- 
driven ions during the source integration steps which starts the neutrals moving after a 
time > T ni has passed for regions near the origin. Once this happens, as we noted above, 
gradients in the neutral density and pressure form, and thermal pressure gradient forces start 
to act on the neutrals during the RG integration step for that fluid. From that point on, 
there is a combination of wave-driving thermal pressure and magnetic forces (through ion 
collisions) on the neutrals during all integration stages of the split-operator method scheme. 
(It is this combination of thermal pressure and magnetic forces acting on the neutrals which 
explains the appearance of both the sound and Alfven speeds in the expression for the neutral 
magnetosound speed Vn ms , eq. [43] .) 



To summarize: the benchmark tests for wave packets show that our split-operator 
method successfully incorporates the relevant physics and accurately follows the evolution 
of both propagating and diffusive flows traveling perpendicular to the magnetic field, over 
many orders of magnitude in both the temporal and spatial scales. 



3.2. Early Evolution of Shocks Caused by a Cloud-Cloud Collision 

We now consider a more dynamic test which shows that the split-operator method is 
also suitable for multifluid shock waves. RC07 considered the collision of two identical clouds 
and found analytic solutions (using Green functions) which describe the early-time behavior 
of the ensuing disturbances in the ion-electron fluid. For the parameters RC07 considered, 
the collision produced forward- and reverse J-shocks in the neutral gas. They referred to the 
disturbances in the ion-electron fluid as "driven waves" because they are driven by frictional 
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coupling to the neutral flow. Their results describe initial stages in the formation of magnetic 
precursors on the J-shocks. 

The flow geometry in RC07 is the same as the geometry in this paper, with fluid velocities 
perpendicular to the magnetic field. The clouds are identical and semi-infinite; the collision 
takes place when their free surfaces coincide at the plane x = at t = 0. Prior to the collision 
the charged and neutral fluids in each cloud move together with t> n (x,0) = V{(x, 0) = v(x), 
where 

, . J —10 km s -1 for x > , , . 

V W = \ + io km s- 1 for x < . { } 

Each cloud has an initial neutral density n n (x, 0) = n n o = 2xl0 4 cm -3 , ion density nio(x, 0) = 
riio = 5.72 x 10 -4 cm -3 , fractional ionization X{(x, 0) = riio/n n Q = 2.86 x 10 -8 , temperatures 
T n (x, 0) = T{(x, 0) = 10 K, and magnetic field B(x, 0) = B — 50 /zG. The ion magnetosound 
speed is Vi ms = Via = 912 km s -1 in the undisturbed charged fluid. 

We adopt the same source terms as RC07, who set S a — G n — A n = and assumed that 
momentum transfer is due solely to elastic ion- neutral scattering. RC07 assumed further 
that ri n is independent of the relative velocity between the ions and neutrals. To make a 
legitimate comparison with the RC07 solution, we drop the term oc \v n — t>i| 2 on the right- 



hand side of equation (14) when calculating F n . The ion-neutral mean collision time in the 
clouds is then r in = 1.26 x 10 -2 yr, and the neutral-ion drag time is r ni = 3.52 x 10 4 yr. 
Because r n ; is much longer than the times considered by RC07, the motion of the neutral 
gas is unaffected by coupling to the ions and magnetic field. 

Figure [5] shows the multifluid MHD split-operator code results for the colliding clouds 
test at t — 0.1r in = 1.26 x 10~ 3 yr. At this very early time (<C T in ), the ions and magnetic 
field have yet to be affected by friction from the neutrals. Two discontinuities in the magnetic 
field and charged fluid, symmetric about the origin, travel outward at the ion magnetosound 
speed V ims , with their fronts located at |xidi SC | = 3.63 x 10 12 cm. There are also two symmetric 
neutral shocks propagating away from x = at a speed i> ns hk — 3.34 km s -1 . The neutral 
shock fronts are located at |x ns hk| = 1-33 x 10 10 cm so they are not visible on the scale of 
Figure [5] Also shown (dash-dot curves) in each panel of the figure are the analytic solutions 
of RC07. The split-operator results agree very well with the RC07 solution, with a relative 
difference in the magnetic field that is < 0.015 behind the discontinuities. 

Results for the colliding clouds test at t — 5r in = 6.29 x 10 -2 yr are presented in 
Figure [6j The shock fronts in the neutral fluid continue to travel away from the origin at 
3.34 km s -1 , and are located at |x ns hk| = 6.63 x 10 11 cm. At this stage of the evolution, 
pronounced precursors in the ion velocity and the magnetic field extend to much greater 
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Fig. 5. — Colliding clouds at time t = 0.1ri n = 1.26 x 10~ 3 yr. Top panel: Velocities of the 
neutrals (solid line) and ions (circles). Also shown is the approximate analytic RC07 solution 
(dash-dot line) for the ion fluid velocity. Bottom panel: Magnetic field increase above the 
initial value in each cloud AB = B(x,t) — B , relative to the initial field strength (circles). 
The dash-dot line shows the RC07 solution for the magnetic field. 
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distances beyond the neutral shock fronts. At the time shown (> r in ), collisional drag from 
the inflowing neutrals is having a pronounced effect on the precursors. Notably, they continue 
to be led by discontinuities in the velocity and magnetic field heading outward at the ion 
magnetosound speed V[ ms , but the jumps at the discontinuities are significantly diminished 
in magnitude and strength compared to the earlier time (Fig. [5]). The RC07 solutions are 
also plotted in Figure [6j Once again the split-operator results agree extremely well with the 
analytic solutions. 

The solutions at time t = 50r- m = 0.629 yr are shown in Figure [7] along with the 
RC07 solutions. The agreement between the split-operator results and the analytic solutions 
is again very good. At this time (>» r in ) the ion and magnetic field discontinuities have 
virtually vanished, with the precursors ahead of the neutral shocks having profiles that are 
smooth and continuous. Inside the precursors the motion of the ions is dictated by a balance 
between the collisional drag from the neutrals flowing toward the origin and the outwardly 
directed magnetic pressure gradient. This is consistent with the RC07 solution, which found 
that on these time scales the solution for the charged fluid motion tended toward a force-free 
diffusion mode. 

Figure [7] also contains insets showing the fluid velocities and magnetic field on an ex- 
panded scale which resolves the region between the neutral J-shocks. The vertical dashed 
lines in the insets mark the location of the forward- and reverse shocks at that time. The 
RC07 solutions agree well with the split-operator numerical simulation, even on these greatly 
magnified scales, with the numerical solution for the magnetic field differing from the RC07 
solution by less than 0.3% within the inset region. This is an especially interesting result: 
to be able to use the method of Fourier transforms to obtain their analytic solutions, RC07 
were forced to assume that the neutral gas density is uniform — even in the region between 
the two neutral shock fronts — despite the fact that for strong (i.e., large Mach-number) 
adiabatic shocks the density behind a shock front increases by a factor of 4 (for 7 = 5/3, 
see, e.g., Zeldovich h Raizer 1966 or Spitzer 1978). RC07 argued on physical grounds that 
their analytic solution would nevertheless be valid. The excellent agreement between the 
analytic solutions and numerical results between the two shock fronts (where the density of 
the neutrals does indeed have a fourfold increase) confirms that their reasoning was correct. 

The colliding clouds test results and their agreement with the analytic solutions of RC07 
show that our operator-splitting scheme faithfully reproduces the physics of dynamic MHD 
multifluid astrophysical flows. The split-operator code successfully models the evolution of 
shocks and discontinuities in both the neutral gas and the charged fluid, having no difficulty 
in handling the interaction between the fluids which leads to the formation of a magnetic 
precursor. 
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Fig. 6. — Colliding model clouds at time t = 5ri n = 6.29 x 10 2 yr. All normalizations and 
symbols have the same meaning as in Fig. [5} 
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Fig. 7. — Colliding model clouds at time t = 50r; n = 0.629 yr. All normalizations and 
symbols have the same meaning as in Fig. [5] The insets show values about the collision 
point at x = 0, including the region between the leftward and rightward traveling neutral 
shock fronts, each located by the vertical dashed lines. 
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3.3. Shocks with Mass Transfer 

As a final test we present a model which allows for the transfer of mass between the 
neutral and charged fluids in shocks. We introduce a nonzero source term for the conversion 
of ion-electron mass to neutral mass, with 

S n = m x (a DR ni 2 - C C R™n) > ( 55 ) 

where ( CR is the cosmic-ray ionization rate and a DH is the rate coefficient for dissociative 
recombination of molecular HCO + with electrons; note that in this expression we have 
used n e — ni, which follows from ([7]). In this test we use the representative ionization rate 
C C r = 5 x 1CT 17 s _1 (Dalgarno 2006) and take a DR from the UMIST Astrochemistry Database 



( pit t p : / /www . udf a . net[ Millar, Farquhar, & Willacy 1997), 

7 /300 K\ - 69 , 1 
a DR = 2.4 x 10~ 7 l—f—J cm s • ( 56 ) 

Mass exchange also adds to the momentum source term, so that 

F n = F a>el + F njine i, (57) 

where 

^n,inel = m (a^n^Vi - CcR^n) ■ (58) 

In this test we again take G Q = A n = 0. Neglecting radiative cooling means that the 
temperatures in this test will be much larger than a realistic model with cooling included. 
However it allows us to isolate the effects of mass-transfer between the neutrals and the 
ions as a test of the split-operator code. The rate coefficient for dissociative recombination 
depends on the electron temperature, which is not calculated in the current version of our 
code. For expediency we set 

T e = max(T n , 0.1571) , (59) 

in rough agreement with realistic simulations of steady multifluid shocks (e.g., see Figs. 1-3 
of Draine et al. 1983). 

Initially this model has a discontinuity at xq = 1.12 x 10 16 cm. For x < x$ the neutral 
gas and ions are uniform and flowing in the +x-direction with a common supersonic velocity, 
while for x > xq the matter is uniform and stationary. The other initial conditions are listed 
in Table [TJ The initial ion density was calculated from the expression 

ni = fCcRM 1/2 ; (60) 
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Fig. 8. — Mass transfer test at t — 632 yr. In each panel the initial state is shown as black 
dashed lines. Also displayed are the locations of a neutral contact discontinuity (vertical line 
with squares) in the MT model and an ion contact discontinuity (vertical line with triangles) 
in the NMT model, (a) Neutral density, (b) Ion density (curve with circles) in the MT 
model and NMT model (curve with diamonds). The quasi-mass-equilibrium relation ( |60| is 
the dash-dot curve, (c) Velocities of the neutrals (solid curve) and ions (curve with circles) 
in the MT model, and the ions in the NMT model (curve with diamonds). 
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Fig. 8. — Cont. (d) Magnetic field in the MT model (curve with circles) and the NMT model 
(curve with diamonds). 



-34- 



which assumes that the creation and destruction rates of ions in eq. (55) are equal. 



Figure [8] presents the mass transfer (MT) test solution at t — 632 yr. For contrast we 
also show results for a model with no mass transfer (NMT) which is identical otherwise. 
At the relatively early time shown in the figure (<C r n i ~ 3 x 10 4 yr) the neutral gas is 
virtually unaffected by ion-electron drag. The effect of mass transfer on the neutral gas is 
also imperceptible in the MT model for the following two reasons: (i) even if all of the ions 
and electrons were to recombine, the added mass to the neutrals would be negligible because 
Pi < 10~ 6 p n ; (ii) the rate of change of the neutral density caused by cosmic-ray ionization is 
orders of magnitude smaller than the density change caused by advection. 

In Figure [SJi it is evident that the collision of the inflowing and stationary material at x 
has resulted in two neutral J-shocks: a left shock located at XL, s hk = 2.60 x 10 16 cm traveling 
with a velocity t>L, s hk = +7.41 km s -1 , plus a right shock at XR jS hk = 3.93 x 10 16 cm with 
velocity t>R iS hk = +14.1 km s _1 . The figure also reveals a contact discontinuity in the neutral 
gas at x n con = 3.22 x 10 cm. The neutral density profile is plotted for the MT model; the 
neutral density in the NMT is indistinguishable from that of the MT model. 

The effect of mass transfer on the ion fluid in the MT model can be seen in Figure 
[8^. Between the shock fronts the neutral gas is heated, with a corresponding increase in 
the temperatures of the ions and electrons. This results in a decrease in the number of 
ion-electron recombinations there (see eq. (56] ). At the same time, shock compression of the 



neutral gas leads to a greater number of cosmic-ray created ions (see eq. 55 ). The interplay 
of these two effects leads to ion densities between the shock fronts which are a factor 5 
greater compared to the NMT model. Not only does the NMT model have a significantly 
smaller ion density in the shocked region, it also has an ion contact discontinuity located at 
^i,con — 3.32 x 10 16 cm. There is no ion contact discontinuity in the MT model because the 
source integration steps in the algorithm completely overwhelm the RG step (Fig. [T]). That 
is, the ion density profile is determined almost entirely by chemistry rather than advection. 

Table 1. Initial Conditions for the Mass Transfer Test 

v n (= Vi) n n rii T n (= 7] = T e ) B V ims 

x < x : 20 km s^ 1 2.5 x 10 4 cm~ 3 8.12 x 10~ 4 cm" 3 15 K 50 pG 765 km s" 1 

x > x : 2 x 10 4 cm' 3 6.31 x 10~ 4 cm" 3 10 K 25 fiG 434 km s" 1 
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There is a small jump in the ion density at the neutral contact discontinuity caused by the 
jump in the neutral density, which is proportional to the rate of ionizations per unit volume. 
Increased Ti and T e (due to significant drift speeds between the ions and the neutrals [see 
eqs. [8] and (59]) inside the magnetic precursors upstream from the shock fronts reduces the 
rate of ion-electron recombinations there. This raises the ion density inside the precursors 
to values noticeably larger than those in the NMT model. 

The width of the magnetic precursors depends on the amount of ion mass loaded onto 
the magnetic field lines there, as a comparison of the MT and NMT models shows (Fig. [8fc). 
The width of a magnetic precursor scales as Lprc ~ Vjms^Tjn / t^shki where f s hk is the shock 
speed; this result can be derived by balancing magnetic pressure and ion-neutral drag in 
the precursor (Draine 1980) or by setting the precursor crossing time t pre = L pre /v shk equal 
to the ion-magnetic field diffusion time ahead of the front t^g — L pTC 2 /2\d (see eq. 



49d 



Because V| ms = V[a in our models, it follows from equation (31) that L pre oc nC 1 . Hence, 
the model with the greater ion density n\ will have the smaller values of L prc . This explains 
why the precursors extend further in the NMT model than in the MT model. 

Figure |8]i reveals a consequence of the dependence of precursor width on ion density: 
the distribution of magnetic flux is noticeably different in the MT and NMT models. Com- 
pression of the magnetic field by a shock is made less abrupt by having a precursor of greater 
width. This accounts for the smaller magnetic field increase at the left shock in the NMT 
model than in the MT model. (Note, however, that magnetic flux is conserved in both 
models: the total magnetic flux contained in the region from 1.5 x 10 16 cm to 4.5 x 10 16 cm 
[= area under the curve] is the same for both.) A by-product of the mass transfer test, 
then, is that we have demonstrated how loading of mass onto magnetic field lines affects the 
structure of a magnetic precursor and the magnetic field profile elsewhere in a shock. These 
effects feed back into the dynamics, because the magnetic pressure gradient is one of the 
main driving forces in the ion-electron momentum transport equation ([5]). 

The fidelity of the mass-transfer test can be quantified by noting that from the left-hand 
side of the ion mass equation Q we can define an ion advection time 

Tiadv = - = 317 — yr 61a 

t>j \10 ib cm/ \ V{ J 

where L is a characteristic length scale. From the right-hand side of the same equation we 
can also define a recombination time scale 

- - « 66 .Q (I >< IQ-Wn/^V- yr _ 



m ; a 7ii 2 V rii J \300 K 
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and a cosmic-ray ionization time 



2 x 10 4 cm- 3 \ /5 x 10~ 17 s" 1 



T10R S ^ = 63 ' 4 (2 x 10* cm- J l." V" ) {" "Cc» " J ^ <61C) 

(see eq. [55]). For the physical conditions in the MT model, r ia d v > T irec , Ticr on length scales 
L icq >5x 10 15 cm. For model ages greater than r ire f and ticr and length scales greater than 
L ieq , ion advection can be ignored, and there should be approximate equality between the 
rates of mass creation and destruction. That is, quasi-mass-equilibrium with S n m should 
occur with n\ given by equation (60). It is therefore expected that the ion density in the 
MT run for the time displayed should be set by the condition of quasi-mass-equilibrium in 
regions having length scales > L ieq . 

To test this hypothesis, the values of ri\ predicted by the quasi-equilibrium equation 



(60) are also plotted in Figure [8p (dash-dot curve). It is seen that the predicted values do 
match the actual MT model results well in most places, except at the shock fronts where 
the assumption of large length scales becomes invalid. Away from the shock fronts the 
agreement between the MT model and equation (60) is generally quite good. For instance in 
the region between the shocks, including the region about the neutral contact discontinuity, 
the relative differences of the quasi-equilibrium and model values for n\ range from 3 x 10~ 3 
to 0.01. The very good agreement between the values predicted by equation (60) for the ion 
density with the actual results of the fully dynamical MT model (in the regions where the 
mass-quasi-equilibrium approximation is valid) illustrates the accuracy of the split-operator 
scheme for multifluid shocks with mass transfer by ionization and recombination. 



4. Summary 

Because many protostellar outflow shocks are much younger than the time scale r n i to 
accelerate the neutral gas, it is likely they are not steady flows. A time-dependent treatment 
is therefore generally required to model outflow shocks and their emission. In this paper 
we have presented a method for modeling perpendicular, time-dependent, multifluid MHD 
shocks using operator splitting. This scheme splits the time integration of the evolutionary 
equations into separate steps, one dealing with the solution of independent homogeneous 
Riemann problems for the neutral and charged fluids using Godunov's method, plus other 
steps where only the equations describing interactions between the fluids are integrated. Our 
method exploits the fact that the thermal pressure of the charged fluid is usually negligible 
compared to its magnetic pressure. Neglecting the thermal pressure reduces the number of 
characteristic waves in the ion-electron fluid to three, the same number as in one-dimensional 
gas dynamics (Toro 2009). We have shown that under these circumstances the MHD Rie- 
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mann problem can be solved exactly. Using this exact solution, we have constructed an 
approximate MHD solver which is used in the Riemann-Godunov integration of the charged 
fluid in our split-operator scheme. The similarity in structure of the Riemann problems for 
both fluids in our split-operator scheme makes it straightforward to adapt well-established 
gas-dynamic numerical techniques (TVD slope limiters, other data reconstruction methods 
such as WENO, etc.) for use in the RG integration of the charged fluid. The symmetry in 
the RG problems for both the neutral and charged fluids also greatly simplifies the numerical 
coding. (In fact, following our scheme it should not be too difficult to modify an already ex- 
isting non-magnetic single-fluid RG hydrodynamics code to one for use in multifluid MHD.) 

Since the inertia of a fluid and the wave modes it supports are fundamental to the 
solution of the Riemann problem, the multifluid split-operator method outlined here has no 
difficulty dealing with flows involving MHD waves or transients in the charged fluid. Several 
tests spanning a wide variety of time- and length scales were performed to demonstrate the 
versatility and accuracy of our algorithm. The numerical results are in very good agreement 
with analytic solutions in all of our benchmark tests. The latter include a model for the 
formation of MHD shocks resulting from the collision of two identical clouds, a problem 
solved analytically by RC07; the split-operator code successfully captures the ion density 
and magnetic field transients which propagate away from the collision surface at early times. 
Our code also correctly reproduces the multifluid shocks with magnetic precursors which 
develop at much later times, when the charged fluid is force free and evolves diffusively. 

Another benchmark test involved MHD shocks with transfer of mass between the neu- 
tral and charged fluids. A numerical model with mass transfer has significantly greater ion 
densities behind the shocks and in the magnetic precursor compared to a model with no mass 
transfer. The enhanced loading of ion mass onto magnetic field lines in the precursors of the 
mass transfer model results in precursors of smaller width, and there is a corresponding dif- 
ference in the distribution of magnetic flux between models with and without mass transfer. 
Analysis of the characteristic time scales for cosmic-ray ionization and ion-electron disso- 
ciative recombination in the mass transfer model suggested that for sufficiently large length 
scales the ion density should be close to the value that would be predicted when ionization 
and recombination balance one another. The ion density profile in the mass transfer model 
is indeed well fit by the quasi-equilibrium relation, except in the vicinity of the shock fronts 
where the criterion of large length scales (or, equivalently, negligible ion mass advection) is 
violated. 

This work was supported by the New York Center for Astrobiology, a member of the 
NASA Astrobiology Institute, under grant #NNA09DA80A. Code development and model 
runs were performed on the facilities of the Computational Center for Nanotechnology In- 
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novations (CCNI), which is partnered with Rensselaer Polytechnic Institute. 



A. The Exact MHD Riemann Solution 



The three-dimensional MHD Riemann problem for an arbitrary flow geometry and a 
plasma with nonzero thermal pressure is probably too complex to be solved exactly (Brio 
& Wu 1988; Torrilhon 2003). Because of this, several approximate MHD Riemann solvers 
have been developed (e.g., Dai & Woodward 1994, 1997; Ryu & Jones 1995; Balsara 1998; 
Li 2005; Mignone 2007). However for the special case studied here — one-dimensional flow, 
perpendicular geometry, and no thermal pressure — the MHD Riemann problem for the 
charged fluid has a simple, exact solution. We give it here. 



A.l. Characteristics 



The Riemann problem for the charged fluid is to solve equation (28), 

, 



(Al) 



dt dx 
subject to Riemann initial conditions, 

U\ = [pf, p\v\, B L ] T ifx<0 

uf ee [pf, p*y\ b r ] t if x > o 

where the constant vectors Uf and U\ are the initial states on the right and left half planes 



C/i(x,0) 



(A2) 



The first step is to rewrite eq. (Al) in the equivalent form 



at ox 

and examine the eigenvalues and eigenvectors of the Jacobian matrix, 



(A3) 



dU u 







-vj 



1 

2Vi 5/47T 



The eigenvalues are 



-Bvi/pi B/pi Vi . 
{A_, A , A + } = {wi-y iA , v h Vi + V iA } 



(A4) 



(A5) 
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where is the ion Alfven speed (eq. [31] ). The corresponding right (column) eigenvectors 
are 

(A6a) 
(A6b) 



R = [pi, Pi {vi - V iA ) , B\ 
Ro = [ph PiVu 0] T , 



and 



R + = [a, Pi (Vi + V iA ) , B] T . 
The solution also depends on the left (row) eigenvectors, which are 

Pi ' P\ B 
1 



2vu 



0. 



Pi 



B 



and 



2^iA 



^i 1 ViA 
Pi' Pi' 5 



(A6c) 

(A7a) 
(A7b) 

(A7c) 



The left- and right eigenvectors are orthonormal. 



Introduce the characteristics curves, 



dt 



A 



x m (t), defined by 
m = -, 0, +. 



(A8) 



For Riemann initial conditions, the characteristics emanating from the initial discontinuity 
at x = are straight lines (Fig. [9]). Each represents a discontinuity in the solution which may 
be a shock wave, rarefaction, or contact discontinuity. Since discontinuities propagate along 
characteristics, and the initial data are uniform to the left and right of the origin, Figure [9] 
shows immediately that Ui — Uf to the left of A_ and Ui = Iff' to the right of A+. The 
problem therefore reduces to finding Ui in the "star region" between A_ and A+. Let p* L , v * L , 
and B* L denote the 3 unknowns in the star region between A_ and Ao and p* R , v* R , and B* R 
be the 3 unknowns between Ao and A + . These are found by (i) classifying each characteristic 
as a shock, rarefaction, or contact discontinuity; and (ii) joining the solutions in such a 
way that the Rankine-Hugoniot (RH) jump conditions are satisfied across shocks and the 
generalized Riemann invariants are conserved across rarefactions and contact discontinuities. 
For a detailed discussion of the underlying theory, see the monograph by Toro (2009). 



A. 2. The An Characteristic 



The classification of Ao depends on the "eigenvalue gradient," 

<9A d\ d\ 



VA 



at/a' dtV dU Q 



(A9) 
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It is easy to show that 



VA 



P, 



[-v i} 1, 0] 



(A10) 



and hence that 

VXo-Ro = 0. (All) 

According to the theory of hyperbolic]^] PDEs, equation (All) establishes that Ao is always a 
contact discontinuity Conservation of the Riemann invariants across the mth characteristic 
implies that 

dU x dU 2 dU 3 



R-ml Rm2 Rm3 

(see Toro 2009, §2.4.4) and writing this out for m = gives 

dpi d (piVi) dB 
Pi piVi 



Integrating equation (A13) across the contact discontinuity yields 



V * L = V * R 



and 



B 



B 



*R 



B*. 



(A12) 



(A13) 



(A14) 



(A15) 



Thus the ion density may undergo a jump across the contact discontinuity but the other fluid 
variables are continuous, in precise analogy with contact discontinuities in gas dynamics. 
Notice that the number of unknowns in the star region has now been reduced to four: v*, 
B*, p* L , and p* R . 



A. 3. The A i and A Characteristics 



The eigenvalue gradients for the other characteristics are 



VA d 

from which it follows that 



Pi 



Pi 



- =F 1 

' Pi' + 



VA± • R± = =f-ViA. 



(A16) 
(A17) 



1 If B ^ the eigenvalues of A are real and nondegenerate. This guarantees that eqs. (A3 1 are strictly 
hyperbolic (e.g., Jeffrey 1976). 
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It seems reasonable to assume on physical grounds that B ^ everywhere for t > (no vac- 



uum) so that the RHS of expression (A17) is always nonzero. Then the theory of hyperbolic 
PDEs establishes that A_ and A + are never contact discontinuities. The A_ wave is a shock 
wave if the pressure in the star region exceeds the pressure in the left region (B* > B h ) 
and a rarefaction wave otherwise. Similarly, A + is a shock if and only if B* > B R . These 
conclusions follow from the fact that shocks are compressive and rarefactions are not, plus 
our assumption that the pressure of the plasma is entirely magnetic. For a given set of initial 
conditions, only one combination of shocks and/or rarefactions at A_ and A + gives a solution 
which satisfies the matching conditions across all three discontinuities. In £ A.4 we give a 
simple criterion for identifying the correct combination. 

The matching conditions for shock waves are the RH jump conditions, which require J-\ 
to be conserved across the shock front in a frame comoving with the shock. We omit details 
of the derivation and simply give the results. For a shock at A_ (a "left shock"), we find 
that the ion velocities on opposite sides of the shock front are related by 



where 



v\ + /ls m 



(P*-P L ) ( B L \ 



1/2 



and P{B) = B 2 /8ir is the magnetic pressure. Similarly, for a right shock we find 



where 



/rs (B* 



(P* - P R ) 



PI 



~B* 



1/2 



(A18) 



(A19) 



(A20) 



(A21) 



The matching conditions for rarefactions are governed by the Riemann invariants. Car- 



(A22) 



rying out steps analogous to the derivation of eq. (A13), we find that 

dpi d (piVi) dB 
Pi pi (vi ± V iA ) B ' 

where the upper sign corresponds to A+. Equating the first and third terms simply gives 
flux freezing in the charged fluid. This implies that the magnetic field and hence Via may 
be viewed as functions of the ion density alone. Knowing this, we equate the first two terms 



in eq. (A22) to obtain a differential equation for v-{. 

dvi = ±V iA (pi) 



dpi 



(A23) 
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Integrating eq. (A23) across A_ gives the matching condition for a left rarefaction: 

< = ^i L + /LR (#*), 

where 



"2^ 



1 



5 1 



*\ 1/2 



(A24) 



(A25) 



and I^a ^ s the i° n Alfven speed in the left region. For a right rarefaction one finds similarly 
that 

v* = v* + f RR (B*), (A26) 

where 

B* 



1/2 



(A27) 



A. 4. Flow Configuration and Solution 

Let RR, SR, RS, and SS denote the four possible flow configurations, where "RR" is 
a flow with two rarefactions, "SR" a left shock plus a right rarefaction, and so on. If the 
configuration was known a priori, one could write the matching conditions across the A_ 
wave, 

vf = v\ + /l (B*) , (A28) 

and the A + wave, 

v* = v* + f R (B*), (A29) 
by choosing f L and f R appropriately from the functions f^s, /lr, etc. Subtracting equation 



(A29) from (A28) gives 



v? + h(B*)-f R (B*) = 0, 



(A30) 



which is a nonlinear equation for B*. Once B* has been determined by solving equation 



(A30), the velocity in the star region follows from equation (A28) or equation (A29). Finally. 



the matching conditions for shocks and rarefactions both require flux freezing in the plasma, 
and this determines the density in both halves of the star region: 



R* 



and 



Pi 



*R 



B* 



(A31) 



(A32) 
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To implement the steps above it remains only to identify the unique flow configuration 
which satisfies the matching conditions across all three characteristics, an exercise in the 
process of elimination. For example, suppose that the flow is of type RR, which requires 
that B* < B L and B* < B R , for the case v\ - vf > (the flows collide). Equation (A30) 
has a solution only if /l — /r < 0; but examination of equations (A25) and (A27) shows 
that /lr — /rr is strictly positive if B* < B L and B* < B R . We conclude then that the 
RR configuration never occurs when vf — v R > 0. Similarly, consider a flow of type SS, 
which has B* > B L and B* > B R , for the situation v\ — v j < (the flows diverge). For 



that situation, according to equation (A30) there can only be a solution if /l — /r > 0. 
Equations (|A19|) and (|A2l|) reveal that, for B* > B L and B* > B R , f L - f R is always < 0. 



We conclude then that the SS configuration cannot occur when vf — v R < 0. Straightforward 
but lengthy analysis of the other possibilities shows that the configuration depends only on 
three dimensionless parameters, 

— y R 

e=^7T7^ ; (A33a) 



9 = B R /B 1 



and 



^ = \lpf/p\ 
The flow classification is given in Table [2j where 



r 

T 



1 



6V2 
1 



(A33b) 
(A33c) 

(A34a) 

(A34b) 

(A34c) 
(A34d) 



Once the flow type is known, which of the expressions (A19), (A21), (A25), and (A27) 
should be used as the functions fh{B*) and f R (B*) in equation (A30) is then also known. 
The value of B* can then be determined to any desired level of accuracy from equation 
(A30) using an iterative numerical technique such as Newton's method (Aktkinson 1989; 
Press et al. 1996). We note that for flows with both left and right rarefaction waves (type 
RR), the ion mass densities in the star region p* L and p* R must be > 0; it follows from 
equations (A31) and (A32) that the magnetic field in the star region B* must also be > 0. 
This non-vacuum (or, positivity) condition imposes a lower limit on the value of £, the 
dimensionless velocity difference between the left and right states of the initial discontinuity 
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(A33a). Using equations (A25), (A27), and (A30), one finds that this condition is satisfied 



so long as £ > £ vac , where 

£ vac = -2^1 + . (A35) 

As demonstrative examples, we consider two representative MHD Riemann test prob- 
lems. In the first problem, the initial state of the plasma and magnetic field has B L = 50 pG, 
B R = 25 /iG, v\ = 200 km s" 1 , t> ; R = 0, and p\ = p R = 2.51 x 10~ 26 g cnr 3 (for a cloud with 
m,\ = 25m p , this corresponds to a number density n\ = 6 x 10~ 4 cm -3 ). For these parameters 
B L > B R , and T = 0.321 > 0. Examination of Table [2] indicates that the flow configuration 
for this first test will be of type RS, a left rarefaction with a right shock. That this is so can 



be seen in Figure 10 which shows the results for this problem. The initial state of the plasma 
and magnetic field are displayed as dashed lines, and the solid curves are the MHD Riemann 
solution at the time t = 0.1 yr. The right shock is located at x s hk = +2.07 x 10 14 cm, and 
the head of the left rarefaction wave is at x rw = —2.49 x 10 14 cm. A contact discontinuity is 
located at x c = +8.23 x 10 13 cm; that is the position the initial discontinuity in the charged 
fluid's magnetic flux-to-mass ratio B j p\ has moved to by the time shown. 



In the second MHD Riemann test problem we have an initial state with B L = 45 pG, 
B R = 50 pG, v\ = —200 km s -1 , v R = +200 km s _1 , and the same constant ion density 
as in the first test. This state has B h < B R , v\ < v R , and $ = —0.385 < 0. According 
to Table [2] this should be an RR flow. This is indeed the case, as can be seen in Figure 11 



which displays the exact MHD Riemann solution for this model at t = 0.1 yr. At the time 
shown, the head of the left rarefaction wave is at XL, rw = —3.14 x 10 14 cm, and the head 
of the right rarefaction wave is at XR irw = +3.42 x 10 14 cm. There is also a charged fluid 
contact discontinuity at x c = —1.70 x 10 13 cm. 



B. An Approximate MHD Riemann Solver 
B.l. The Star Region 

The Riemann solution found in Appendix [A] can also be used to numerically solve the 



MHD flow problem of equation (28) using Godunov's method (Godunov 1959; Toro 2009). 
However, doing so can be computationally expensive, because it uses an iterative method to 



solve for B* from equation ( A30 ) . 
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Table 2: 






Classification 


of Solutions for the Charged Fluid 








Initial Conditions 


Flow Type 


Note 


v t 


- vf 


> : 


B L 


< B R , * < 


SS 


left shock + right shock 








B L 


< B R , ^ > 


SR 


left shock + right rarefaction 








B L 


> B R , r < 


SS 


left shock + right shock 








B L 


> B R , r > 


RS 


left rarefaction + right shock 


v t 


- vf 


= : 


B L 


< B R , * > 


SR 


left shock + right rarefaction 








B L 


> b r , r > o 


RS 


left rarefaction + right shock 


v\ 


- vf 


< : 


B L 


< B R , $ < 


RRt 


left rarefaction + right rarefaction 








B L 


< B R , $ > 


SR 


left shock + right rarefaction 








B L 


> T < 


RRt 


left rarefaction + right rarefaction 








B L 


> T > 


RS 


left rarefaction + right shock 



tThe RR-state can exist only if -2(1 + 6/5)V& < v\ - vf < (see eq. A35 ). 
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It is desirable then, to instead use a more efficient approximate solver to calculate the 
variables in the star region (B*, v*, p* R , and p* L ) between the characteristics A_ and A + (see 
Fig. [9|. This can be done by assuming that the flow motion is approximately linear, that 
is, a wave, and then using the characteristic relations derived from the Riemann invariants 
(eq. |A22| ) to connect the variables in the star region to those in the adjacent L or R region 
of the Riemann problem (see Fig. Km. But this is just what was done in deriving equations 



(A25) and (A27). Solving those two equations for B* and v* in terms of the L and R region 
variables yields 



B* 



where 



c\ + Cf 

C\vf + Cfv\ + 2C\Cf (yfW - ^B^j 



C\ + Cf 



C\ 



V L 
s/B 1 



c R 



With these quantities now known, p* L and p* R can be calculated directly from equations 



'B 



R 



(Bl) 
(B2) 

(B3) 



(A31) and (A32), respectively. 



B.2. Calculation of Fluxes at the Discontinuity 

To solve the MHD Riemann problem for a state having the discontinuous initial con- 



ditions (A2) using the method of Godunov (1959), we need to calculate the the flux vector 
^(0) at the interface separating the two states, x = (or, equivalently, along the ray with 
the similarity variable x/t = 0, which corresponds to the t-axis in Fig. [9l. To calculate 



^i(O) from equation (20), we need the array of variables C7j(0) (eq. IT 



A very thorough discussion of how to obtain the state variables at the interface x = 
for the one-dimensional gas dynamic Riemann problem using Godunov' s method is pre- 
sented in Ch. 6 of Toro (2009); the fundamental ideas described there - - in which the 
Riemann- Godunov flux is determined by examining the characteristic wave structure at the 
discontinuity — remain valid and carry over directly to the MHD Riemann problem that we 
are concerned with here. The first step in the assignment of Ui(0) is made by finding the 
values of the velocity v* (= Ao) of the contact discontinuity and the magnetic field of the 



star region B* from equations (Bl) and (B2) 



For v* > 0, the contact discontinuity is moving to the right. If B* > B L the left wave 
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is a shock wave with velocity 

s f = v f - Q\/ p \ (B4) 

in the stationary frame. This relation was obtained by applying the RH conditions in a 
frame moving with the shock front. The quantity 

(W -(W (B5) 

V* — V ; 

is the ion mass flux across the front. If s\ > 0, the shock is moving to the right and 
Ui(0) = U\. If s\ < 0, the shock is instead heading to the left, and 1^(0) = U* L = 
[p* L , pfv* ,B*] T 

If B* < B L , the left wave will instead be a rarefaction. As is widely-known (e.g., 
Courant & Friedrichs 1948; Landau & Lifshitz 1959; Zeldovich & Raizer 1966; Toro 2009) a 
rarefaction emanating from a discontinuity has a fan-like phase-space structure (a centered 
simple wave) with a head and a tail. The velocities of the head (hf) and tail (tf) of the left 
rarefaction fan are 

hf = vl-V^ A} tl = v*-V*jt, (B6) 

respectively, where 

V* A L = B * . (B7) 

For h\ > 0, t\ > 0, the entire rarefaction fan is heading to the right and EA(0) = U\ , while 
for h\ < 0, t\ < the entire fan is traveling left and EA(0) = U* L . If h\ < and t\ > 
values inside the rarefaction fan are needed (e.g., see Ch. 4 of Toro 2009). In that situ ation 



Ui(0) = U\ i&n = [pP an , pL fan Wl Lfan , B Uan ] T . Using the Riemann invariant relation (A22) 



to integrate across the A_ characteristic to connect quantities in the left rarefaction fan to 



those in the left state, along with the relations (A5), (A8), and (A31), and setting x/t = 
in the resulting expressions yields 



,,Liai) _ ^ / L 1 t\ R Lfan _ dL f ^ i ^ „Lfan _ Pi ^ / p< 



If f* < 0, the contact discontinuity is traveling to the left. For B* > B R the right wave 
is a shock wave. The velocity of the shock in the stationary frame is 

= vf- + Qf/pf- , (B9) 

where 

( W -(f)V8* 
v* — V ; 
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is the ion mass flux across the shock front. If sf* > the right shock is traveling to the 



right and 17^0) = U 
Ui(0) = Uf instead. 



*R 



[pf 



p* K v* , B*] ; if sf < the shock is heading left, and 



For B* < B R the right wave will be a rarefaction. The right rarefaction fan has head 
(/if) and tail (tf 1 ) velocities 



hf = + V$ , tf = v* + V{f 



where 



v iA — 



B* 



(bh; 



(B12) 



U* R , while, for 
U?. If 



V^pt R ' 

For > 0, tf > 0, the entire fan is rightward traveling and Ui(0) 
hf- < 0, tf- < the complete fan is heading to the left, and therefore Ui(0) 
hf- < and tf 1 > values from inside the right rarefaction fan are requir ed an d t/i(0) 



U 



Rfan 



Rfan 



Pi 



Rfan„,Rfan 



V; 



5 Rfan ] T . Integrating the Riemann invariant (A22J) across 



A. 



to link quantities in the tail of the right rarefaction fan to quantities in the right state, also 
using equations (A5), (A8), and (A32), and taking x/t = in the resulting relations, gives 



Rfan 



,R 



B 



Rfan 



B l 



Pi 



Rfan 



R RRfan 



P?B 



B 



R 



(B13) 



As we have stated earlier, approximate Riemann solvers for the one-dimensional neutral 
gas dynamic Riemann-Godunov problem have a very similar structure to that of the MHD 
solver we present here. In several instances the algorithm for the MHD and neutral gas 
Riemann solvers are basically the same when a suitable substitution is made. For example, 
the relations we derived for quantities in the left and right rarefaction fans (eqs. |B8| and 
B13| ) are identical to the left and right fan relations at x/t = for the density, velocity and 



thermal pressure of an ideal non-magnetic gas given by equations [4.56] and [4.63] of Toro 
(2009) if one makes the substitution P = B 2 /8tt, replaces the speeds of sound in the left and 
right states with and V^, respectively, and uses the fact that the adiabatic index for a 
flux-frozen plasma 7 = dlnP/dlnpi = dln(B 2 /8ir)/dhipi = 2 in Toro's expressions. 

Although a linear approximation was made in its derivation, we find the MHD Riemann 
solver presented here to be quite accurate and robust, even when dealing with shocks. It 
matches the exact Riemann solution for a variety of test problems, including conditions 
likely to be encountered when studying high-velocity shocks and flows in interstellar clouds. 



Figures 10 and 11 show the results of RG calculations (circles) for the two example MHD 



Riemann problems described in § |A.4| The simulations using the approximate MHD solver 



are seen to be in very good agreement with the exact Riemann solutions for those model 
tests. 
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C. Analytic Benchmark Solutions 



We present here a derivation of the analytical small-amplitude wave solutions used to 
test the split-operator models discussed in §§ |3.1.1 - 3.1.3 For the governing fluid equations 
we linearize the system in the standard way, assuming for any physical quantity / that 



f(x,t) = f + Sf(x,t); 



(CI) 



the zero-order value /o is taken to be constant and uniform, and the first-order perturbation 
\8f(x,t)\ <C fo- Inserting these perturbations into equations Q - ^ and retaining terms 
up to first order give the linearized system of equations: 



0_ 

Of 
d 



PnO^ On) 

J_d_ 
T^dt 
d_ 

di 



(fa) 



Of 



(SB) 



(7-1)— | (W 

PnO Ot 

-P i0 ^ (fa) , 
_Bod_ 
A-k dx 

-B - 
dx 



PnOr . PnOr 
0V n H 0V[ 



(SB) 
(S Vi ) 



PiO r PiO r 

Sv n Of; 



(C2) 
(C3) 
(C4) 

(C5) 
(C6) 
(C7) 



where we have used the fact that the charged and neutral fluids are at rest in the zero-order 
reference state (v n0 = t> i0 = 0), and also that F njine i = S n = Si = G n = A n = for all of the 



model tests in § 3.1 Equation (13) was also used to substitute for F n in the perturbed force 



equations (C3) and (C6) 



The internal energy for an ideal gas (11) having a ratio of specific heats 7 was used in 



the energy equation of the neutral gas (j3]) to derive (C4); integrating that latter equation 
and inserting the result into the linearized ideal gas law (|9]), 



SP n Sp n ST, 



P, 



no 



+ q 

PnO 1 n0 



(a 



yields the adiabatic relation between density and pressure perturbations in the neutral fluid, 

7^n 



SP n 



PnO 



-Sp n = C n Sp n , 



(C9) 



where C n is the adiabatic speed of sound (44). 
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Combining equations (C5) and (C7) and integrating we immediately obtain the result 

8 Pi /5B = Pi0 /B , (CIO) 
which simply expresses freezing of magnetic flux in the charged fluid. Inserting this result 



in equation (CI) and using n- x = Pi/rrii, equations (48b), (49b), and (50b) follow directly 



C.l. Solution for Lq <C Lj ms 



In this limit the neutral fluid is unable to respond to the very short time and length scales 
of the perturbations; therefore, for this situation the neutrals are effectively a fixed stationary 
background with 8v a = 5p Q = 0, and the only equations relevant to these particular modes 



are (C6) and (C7). 



The Fourier transform / of any function / is given by 

1 



/(M) 



v/27T J-oo 

where k is the wavenumber. The inverse transform is 

1 



f(x,t)e~ lkx dx 



f(x,t) 



f{k,t)e ikx dk 



(Cll) 



(C12) 



Fourier transforming equations (|C6|) and (|C7|) gives a system of linear ordinary differ- 

(C13) 



ential equations with constant coefficients: 



dt 



Mil/; 



where the vector of transformed variables yjjz^t) = [5vi(k,t),SB(k,t)] T and 



-l/r in -iB k/4iTp i0 
-iB k 



The ODE system (C13) has the solution 



(C14) 



(C15a) 



(e.g., see ch. 3 of Braun 1983), where the eigenvalues and eigenvectors of are, respectively, 

21 V2 



l 



Tu 



k 



(C15b) 
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and 



Si 



± 



h-i- 



B k 
9i± 



l T 



(C15c) 



In these expressions we have used equation (31 ) and the fact that V; ms = Vi\ for all conditions 



of interest in our model interstellar clouds. We also introduced the ion magnetosound cutoff 
wave number 

fc™ = = ; (C16) 



2Vj ms T; n -^ims 

modes with k > k ims propagate as ion magnetosound waves. 



The constants a- 1+ and ai_ in equation ( |C15a ) are determined by applying the initial con- 
ditions 8v{(k, 0) = (ions initially at rest) and the Fourier transform of the initial Gaussian 
pulse in the magnetic field (see eq. [46]), 



6B(k,0) = — — exp 



(C17) 



Doing that yields 
6B{k,t) = B ° ApLG 




exp 



t 



2ri 

i k\ms I k 



+ 



1 + 



[i-(fc ims A) 2 ] 1/2 . 

ikifns/k 

[i - (WW 72 



exp 



exp 



W* 2 



+ iVimskt 



1 - 



V* 2 



1 Vim.q/c^ 




(C18) 



and 



8vi(k, t) 



2y/2[l- (hms/k) 



21 V2 



exp 




exp 



exp 



+ iV[ ms kt 



l/2> 



La 2 k 2 



iVimskt 



k u 



1/2^ 



. (C19) 



5B(x, t) and 8v{(x, t) are obtained by inverse transforming (C18) and (C19). Performing 
the integration in the inverse transforms is greatly aided by noting that because the initial 
perturbation in the magnetic field excites wavelengths primarily in the region ~ Lq, the 
condition Lq <C L ims implies that the contributions in the Gaussian packet are principally 
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from waves with k 3> k ims . Taking the limit (ki ms /k) 2 — > will therefore introduce only a 
negligible error when inverting the transforms. Applying this limit gives the final solutions 



5B(x,t) 



5vi(x, t) 



A 

B —exp 

T / 

Vims-^exp 



t 



2r h 



exp 



T 2 



exp 



exp 



2r, r 



erf 



x -\- V\ ms t 



2r ir 



exp 



— ix 



L G 

2* 



erf 



exp 



2\ 1 



L 



G 



(C20) 
(C21) 



The solution corresponds to left- and right-traveling wave pulses propagating at the ion 
magnetosound speed VJ ms through a background of stationary neutrals. The decay time 
scale for the pulses is 2r- in , a value that is consistent with the decay of individual Fourier 
wave modes at wavelengths <C L ims (e.g., see § 3.1.1 of Ciolek et al. 2004, or § 3.2.1 of 
Mouschovias et al. 2011). Equations ( |48a ) and (48c) follow immediately from inserting 
([C20l and ((0211 into flCTL 



C.2. Solution for L ims <C L G <C L r 



In this limit the neutrals still remain motionless. Also, the collisional drag of the neutrals 
on the ions essentially balances the driving magnetic pressure gradient in the charged fluid 
force equation. The solution in this wavelength regime can be derived from the Fourier- 



transformed magnetic field and velocity relations (C18) and (C19), by rewriting them as 
BqA p L g 



SB(k,t) 



x 



2V2 
1 



exp 



t 



2r ir 
1 



[1 - (^Aims) 2 ] 1 / 2 



+ 



1 + 



[1 - (fc/fc ims ) 2 ]V2 



exp 



exp 



L G 2 k 2 



k 



L G 2 k 2 



l/2> 



(C22) 



and 



5vi(k, t) 



-iV iras A p L G k 



exp 



2V2k ims [l-(k/k ims y] 1/2 V 2T; 



L G 2 k 2 

eXp I ; Vjms/Cjmst 



4 



1 - 



kh 



Lr 2 k 2 



— exp 



I - ^ims^ims^ 



1/2' 



(C23) 
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The condition L G ^> L ims for the modes considered here is equivalent to taking the limit 
k/ki ms — > 0. Expanding (C22) and (C23) in the ratio k/ki ms and keeping terms only up to 
second-order will then suffice to give the principal solution when inverting the transforms. 
Doing this, and using the relations (C16) and (49d) to replace the cutoff wavenumber k[ ms 
with the diffusion coefficient V ad , gives 

B A P 



5B(x,t) 



exp 



1 - 



8vi(x, t) 



(l+4V ad t/L G 2 ) 1/2 

2B A p V ad T in/ Lq 2 

(l + 4P ad r in /L G 2 ) 3/2 

2B A p V ad T in /L G 2 
[l-AV ad r hl /L G 2 f /2 

2D ad ApX 

L G 2 (l + AV ad t/L G 2 ) 
L G 2 (l-AV ad t/L G 2 ) 



-x 2 /L G 2 



1 + 4V ad t/L G < 



2x 2 /U 



1 + 4V ad t/L G 2 

2x 2 /L G 2 
-4V ad t/L G 2 

-x 2 /L 



exp 



cxp 



-x 2 /L 



G 



1 + 4V ad t/L G 2 
-t 



cxp 



x 2 /L 



G 



AV ad t/L G 



(C24) 



2\3/2 



exp 



1 + AV ad t/L G 
f 



2\3/2 



cxp 



exp 



-x 2 /L 



G 



1 - 4V ad t/L G < 



(C25) 



The solution here is one of ambipolar diffusion, in which the charged fluid and magnetic 
field diffuse through a sea of fixed neutrals, as reflected by the first terms after the equalities 
in eqs. (C24) and (C25). The diffusion coefficient T> ad reflects the balance of the magnetic 
field pressure (since V- irn 2 cx B 2 ) and collisional forces (represented by r in ). From expressions 
(C24) and (C25) the characteristic decay time scale r dec for this solution is found to be the 
ambipolar diffusion time r ad = L G 2 /AV ad . 



C.3. Solution for L G > L x 



In this limit, which has length scales ^> L ims (since L nms ^> L ims ) for our clouds, the 
charged fluid continues to move in a force-free manner, For this situation we can make the 
approximation piod(5v{)/dt ~ in the linearized ion momentum equation (C6). Hence, the 
ions will be traveling at the terminal drift velocity 



5v\ 



B T in d 
47rp i0 dx 



(SB). 



(C26) 



Inserting the expression (C26) into the linearized force equation for the neutral fluid (C3) 



and also the linearized magnetic induction equation (C7), then Fourier transforming those 



equations along with the perturbed mass continuity equation (C2) yields the approximate 
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system of ODEs for the wave modes in this wavelength region, 



dt 



where y n = Sp n (k,t), Sv n (k,i) , 5B(k,t) 
3x3 array of constant coefficients 



M n y n 



(C27) 



is the vector of transformed variables, and the 



M, 



—ip n0 k 

-ikC Q 2 /p n0 -ikB /47rp n0 
-iB Q k -V aA k 2 



(C28) 



In typical clouds with (C n /V nms ) 2 <C 1 (see eqs. 43 -|45|), the system (C27) has the 
solution 

y n = a n+ S Il+ e^+ t + a n _£ n _e^-* + a^S^e^ 1 , 



where 



9n± 



9n0 



-V ad k 2 ± iK m i 



1/2 



V ad k 2 = -V th k 2 



(C29a) 

(C29b) 
(C29c) 



are the eigenvalues of M n , and 

9n± 9n±B 



£n0 



1 , i 



Pnok ' (g n ± + ^ad^PnO 
9n0 9iioBq 



1 ,T~ 



V u - 



B 



1 i 

Pnok ' (fl-nO + ^ad^ 2 )PnO. 



are their associated eigenvectors. In the above, 

2tt 2K, 



1 ,-i 



2K 



PnO PnO 



(C29d) 



PnO 



Vrirr 



So 

PnO 



(C29e) 



^rnns VnA 7"ni ^ad 



(C30) 



is the neutral magnetosound cutoff wave number (eqs. |49d , 31 , 45 and 15 were used 
to derive the latter equalities). For k < k nms the modes with eigenvalues g n ± are neutral 
magnetosound waves; these waves decay by ambipolar diffusion of the charged fluid and 
magnetic field with respect to the neutrals. The g n o mode is a neutral pressure-driven 
diffusion mode (e.g., see Figs, lb and Id of Ciolek et al. 2004; or § 3.2.1 of Mouschovias et 
al. 2011), in which the thermal-pressure gradient force of the neutral fluid drives the neutrals 
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through a stationary background of the charged fluid and magnetic field. The balance of 
thermal-pressure and neutral-ion collisional drag in this mode has the neutrals moving at a 
terminal drift speed, resulting in neutral diffusion with the coefficient X> t h defined in equation 
(53). Simplifications resulting in the latter equalities of (C29d) and (C29e) were made by 
noting that for (C n /V nms ) 2 <C 1, \g n o\ <C |<?n±|- The limit k <C k nms was also used. 



The coefficients a n+ , a n _, and a n0 in (C29a) are found by applying the initial conditions 
5p n (k,0) = 5v n (k,0) = 0, and 5B(k,0) given by (C17). Doing this, we have the solutions 
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Inverting the transforms give the solutions for 5p n (x,t), 5v n (x,), and 5B(x,t). The 
inversion is aided by the fact that for the initial perturbation we consider here, with L G ^> 
L nms , the overwhelming bulk of the wave modes contributing to the wave packet have k <C 
k nms . Hence, taking the limit (k/k nms ) 2 — > will introduce only a minor error when inverting 



the Fourier transforms (C31), (C32), and (C33). Doing this yields: 
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where the two functions Q\ and Q 2 are respectively defined by equations (51) and (52). 
The solutions clearly show that for initial perturbation with L G ^> L nms left- and right- 
propagating neutral pulses traveling at the magnetosound speed Vn ms . There is some damp- 
ing iin the pulses because of ambipolar diffusion, occurring on a time scale r dec = 2r ad = 
L G 2 /2V ad . There is also a part to the above solutions that indicate an effect due to pressure- 
driven diffusion of the neutrals, as can be seen by the terms in the above equations involving 
X\h- Those terms show that pressure and density gradients are created in the neutrals (by 
collisional drag from the ions) that can also affect the motion of the neutral fluid. How- 
ever, these effects are small compared to that which results from the magnetically-driven 
collisional drag of the ions on the neutrals. 
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Finally, equations (50a)-(50d) are obtained by combining (CI), (C34) - (C36), and 
(CIO), and using the relations n n = p n /m n and n- x = Pi/rrii. Equation (50e) follows from 
taking the partial derivative d/dx of (C34) and then inserting the result into (C26). 



D. Numerical Convergence 



Here we provide quantitative data verifying the scaling and convergence of the numerical 
algorithm described in this paper. We do this by modeling the evolution of small-amplitude 
disturbances and comparing the numerical model results against exact solutions of the lin- 



earized system of equations (C2) - (C7). Inserting (C8) and (C9) into those equations, and 



then Fourier transforming them in space (eq. |C11 ) gives 
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where y tot = [Sp n (k,t), 5v n (k,t), 5pi(k,t), 5vi(k,t), 5B(k,t)] T is the vector of transformed 
variables for the full system, and 
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This particular ODE system has the solution 
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(Braun 1983, Ch. 3), where g m is the mth eigenvalue of M tot and £ m is its correspond- 
ing eigenvector. The constants a m are determined by the initial conditions. Making the 
subsitution 



9r, 
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is a complex frequency (with real and imaginary parts u> mt R and co m j, respectively) and 
inverting the transform (eq. |C12| ), yields 
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Vtot = [$Pn(x,t), 5v a (x,t), 5pi(x,t) , 5vi(x,t), 6B(x,t)] T is the total vector of perturbed 
variables. 

A single mode with m = I at wavelength A p and wavenumber k p = 2tt/\ p is excited 
from the initial complex perturbation 

y tot (x,0) = b P £ l e^ x+ ^ ) . (D6) 

The constant b p is taken to be real, and Xp is the phase constant of the initial state. Setting 
t = in (D5) and equating it to (D6) reveals that for the former relation to describe the 
evolution of the latter, we must have a m = y/2Trb p 5i rn 5(k — k p )e % ^ . Hence, 



y tot (x,t) = b p £ l e l ^ x ^' t+ ^ . 
£i is also complex, and can be written as 

£\ = £lr + i£u 



(D7) 



(D8) 



where £i jR and £ij are both real vectors. Inserting (D4b) and (D8) into (D7), it follows that 
the real part of the singly-excited eigenmode of the system is given by 



U[y tot (x, t)) = b p [£ iiR cos (k p x - u l>R t + x P ) - £i,i sin (k p x - ui )R t + x P )) e 



rt 



(D9) 



To test convergence we excite a single mode of the multifluid system and compare the 
resulting evolution of our numerical code to the exact solution (D9). A A p is selected and 
the eigenvalues and eigenvectors of the array M tot at that wavelength are calculated using 
widely available numerical routines (e.g., EISPACK, LAPACK) or software packages (e.g., 
MATLAB, Maple, IDL). To ensure that an excited mode always remains linear, the value of 
the numerical constant b p of the initial state (the real part of eq. D6 , or equivalently, eq. 



D9 at t — 0) is chosen so that its magnetic field perturbation has a relative amplitude that 



is 10 that of the background magnetic field. For our convergence test runs we assume the 
same background reference state and conditions for all of the physical variables (n n0 , n i0 , B , 
T n , etc.) as in the models described in § 3.1 All of the speeds (V[ ms , Vn ms , and C n ), collision 
times (r in and r ni ), and characteristic length scales (L ims and L nms ) are therefore also the 
same as for those models. Additionally, we set the eigenmode wavelengths A p equal to the 
widths L G of the Gaussian test packets of § 3^; hence, the underlying physics of the waves 
discussed for those models also applies to the models presented here. 

Our first excited eigenmode model has A p = 9.35 x 10 11 cm, which is less than the upper 
ion magnetosound wave cutoff length L ims . The selected mode is an ion magnetosound 
wave traveling in the +a;-direction with velocity +V;ms- At this wavelength the mode has 
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Ui ; r = 5.84 x 10 4 s 1 and uij = —1.26 x 10 6 s 1 , the latter value corresponding to a 



decay rate r dec = 0.0251 yr = 2r in . Figure [12p shows the exact solution for the ion velocity 
(solid curve) at t = 1.37 x 10~ 4 yr, by which time the wave has advanced by more than | 
of a wavelength from its initial state (dashed curve). Also shown are the results of three 
models having different spatial resolution: the first model (displayed as plus signs) has 16 
mesh cells within each wavelength of the initial perturbation, the second (squares) 64 cells 
per wavelength, and the third (circles) 128 cells per wavelength. For this perturbation the 



numerical time step is determined by the CFL time step (29); the CFL number v = 0.8 in 
each of these models. Figure \TTf ) displays a zoom-in of the same data about x — 0. The 
local relative error in the ion velocity \[vi }exajc t{x,t) — fi, n um(^, t)] /^exactO^ t)\ is presented in 



Figure 12 showing how the numerical accuracy increases with finer mesh spacing. Finally, 



Figure 12 i shows the error measure in the ion velocity across the entire compuational domain 
of N cells 

1 N 

Li(vi) = — ^2 \v i)eKac t{xj,t) - Vi }Xmm (xj,t)\ (D10) 
j'=i 

(curve with triangles) for several different model runs, all with v = 0.8, as a function of 
each model's numerical resolution. The range is from 4 to 512 cells per wavelength of the 
perturbation. Also shown (dash-dot curve) is a reference curve with a logarithmic slope 
= —2, the value expected for a code that is second-order accurate in space and time. Our 
numerical results are seen to be consistent with second-order accuracy. 

The second excited eigenmode is at A p = 3.74 x 10 15 cm, which is in the range L; ms < 
A p < L nms . The ambipolar diffusion mode that exists within this range is chosen for this 
test: at this wavelength that particular mode has u^r = 1.30 x 10~ 26 s _1 and u\j = —8.47 x 



10~ 9 s _1 . Its decay time is therefore r dec = 3.74 yr. Figure 13 1 displays the initial state and 
the exact solution for the ion velocity at t — 3.17 yr. Also shown are the results of numerical 
runs having resolutions of 16, 64, and 128 cells per wavelength of the perturbation. For 
these model runs, the upper limit to the maximum stable numerical time is limited by the 
time scales occurring in the source terms (here, the ion-neutral collision time T m ) instead 
of that which would be calculated with a traditional CFL time step with v ~ 0.1 — 1 (see 



2.3). Stability of the time integration for these models thus requires that their time step 
At = Atsj < Tin- In analogy to convergence tests in which the CFL number v oc At/ Ax is 
kept fixed, the models for this test have At/ Ax = constant, such that 

At = (Ay^) Tm ' (D11) 

where iVp (> 4) is the number of cells within one wavelength of the perturbation. A close-up 



showing the results for this eigenmode about the origin is presented in Figure 13 5, and the 
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local relative error is displayed in Figure 13 z. The L\{v\) error measure of this mode can be 
seen in Figure 13 i. Examination of that panel indicates that our code tends to scale slightly 
better than second-order accuracy in this instance. 

The final excited eigenmode has A p = 2.24 x 10 17 cm > L nms . For this test we chose the 
mode that corresponds to a neutral magnetosound wave propagating in the +x-direction. At 
this wavelength u^r = 1.02 x 10 -12 s _1 and u)\j = —7.93 x 10~ 13 s _1 ; hence, the decay time 
by ambipolar diffusion for this mode is Td ec = 3.99 x 10 4 yr. Figure 14 z shows the eigenmode 
solution for the ion velocity at t — 5.03 x 10 4 yr, by which time the wave has moved about j of 
a wavelength from its initial state. Three model runs with different numerical resolution are 
also shown: one with 64 cells per wavelength (squares), another with 128 cells per wavelength 
(circles), and the last with 256 cells per wavelength (diamonds). For these eigenmode runs 
it is again the case that the maximum stable numerical time step cannot exceed T m , so the 
relation (Dll) was used to fix the ratio At/ Ax to the same value in each model to test 
convergence in the usual way. An expanded view of the ion velocity data near x = is 
shown in Figure 14 j, and the local relative error in v\ as a function of position for the three 
models is found in Figure 14:. The error measure L\{v\) as a function of numerical resolution 
is provided in Figure 14 i; the data there show for the large part that our code again exhibits 
second-order accuracy. 
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Fig. 10. — First MHD Riemann demonstration problem. Dashed lines indicate the initial 
state of the charged fluid and magnetic field. The solid line is the exact MHD Riemann 
solution at t — 0.1 yr. Circles are the solution calculated using the approximate MHD 
Riemann solver (some data points have been omitted for clarity). Top panel: magnetic field. 
Middle: ion number density. Bottom: ion velocity. 



-67- 




Fig. 11. — Same as in Fig. 10 but for the second MHD Riemann demonstration problem. 
Data shown at t — 0.1 yr. 
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Fig. 12. — Rightward-traveling ion magnetosound wave with A p = 9.35 x 10 11 cm < L ims at 
time t = 1.37 x 10 -4 yr. (a) Ion velocity. Solid line is the exact eigensolution (eq. |D9| )at 
that time, dashed-line is the initial state. Also shown are numerical code results for three 
models having different numbers of mesh cells per wavelength of the initial perturbation: 16 
cells (plus signs), 64 cells (squares), and 128 cells (circles), (b) Zoom-in of the ion velocity 
data about the origin, (c) Relative error in the ion velocity for the same three numerical 
models, (d) Numerically-determined error measure of the ion velocity Li(vi) as a function 
of the number of cells per wavelength in a model (curve with triangles). The dash-dot line 
is a reference curve with a logarithmic slope of -2. 
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Fig. 13. — Ambipolar diffusion mode with A p = 3.74 x 10 15 cm (L ims < A E 



t = 3.17 yr. All symbols and curves have the same meaning as in Fig. 12 
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Fig. 14. — Right-traveling neutral magnetosound wave with A p = 2.24 x 10 17 cm > L nms 
at t — 5.03 x 10 4 yr. (a) Ion velocity: exact eigensolution (solid curve) and initial state 
(dashed). Also displayed are results for models having different numerical resolutions: 64 
cells per initial perturbation wavelength (squares), 128 cells per wavelength (circles), and 256 
cells per wavelength (diamonds), (b) Close-up of the ion velocity data about the origin, (c) 
Relative error in the ion velocity for the three numerical models, (d) Error measure Li(t>i) as 
a function of model resolution. The dash-dot reference curve has a logarithmic slope = —2. 



